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Abstract 

Additive Runge-Kutta (ARK) methods are investigated for application to the spatially discretized 
one- dimensional convection-diffusion-reaction (CDR) equations. First, accuracy, stability, conservation, 
and dense-output are considered for the general case when N different Runge-Kutta methods are 
grouped into a single composite method. Then, implicit-explicit, N = 2, additive Runge-Kutta (ARK 2 ) 
methods from third- to fifth-order are presented that allow for integration of stiff terms by an L-stable, 
stiffly-accurate explicit, singly diagonally implicit Runge-Kutta (ESDIRK) method while the nonstiff 
terms are integrated with a traditional explicit Runge-Kutta method (ERIv). Coupling error terms are 
of equal order to those of the elemental methods. Derived ARK 2 methods have vanishing stability 
functions for very large values of the stiff scaled eigenvalue, — — 00 , and retain high stability 
efficiency in the absence of stiffness, 2 ^ — > 0. Extrapolation-type stage-value predictors are provided 
based on dense-output formulae. Optimized methods minimize both leading order ARK 2 error terms 
and Butcher coefficient magnitudes as well as maximize conservation properties. Numerical tests of the 
new schemes on a CDR. problem show negligible stiffness leakage and near classical order convergence 
rates. However, tests on three simple singular-perturbation problems reveal generally predictable order 
reduction. Error control is best managed with a PID-controller. While results for the fifth-order 
method are disappointing, both the new third- and fourth-order methods are at least as efficient as 
existing ARK 2 methods while offering error control and stage- value predictors. 
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1. Introduction 


It is oftentimes useful to consider the compressible Navier-Stokes equations (NSE) as evolution 
equations with several driving forces, each having somewhat different characteristics. Typically, one 
distinguishes among terms such as convection, diffusion, and reaction. As such, one often considers the 
more tractable convection-diffusion-reaction (CDR) equations as a prologue to the full compressible 
NSE. 34 In the search for ever more efficient integrators, it is intuitively appealing to seek individual 
integration methods that are ideally suited for specific parts of the governing equations. The individual 
methods are then rolled into a single composite method that, ideally, would be more efficient than any 
individual method applied to the full computation. To accomplish this, one may consider partitioned 
methods. Schemes constructed to take advantage of termwise partitioning of the CDR for integra- 
tion purposes may be called additive methods. 48 Partitioning of the discretized form of the equations 
may also be performed on an equationwise or pointwise basis. 24 There are many different partitioning 
strategies.' 0 Runge-Kutta methods, with their extensive theoretical foundation, allow for straightfor- 
ward design and construction of stable, high-order, partitioned methods composed of arbitrary numbers 
of elemental Runge-Kutta schemes. In addition, they also allow for the direct control of partitioning 
(splitting or coupling) errors. 40 Direct numerical simulation (DNS) and large-eddy simulation (LES) of 
fluid phenomena, with their relatively strict error tolerances, are prime candidates for such methods. 
The need for these strategies is by no means limited to Navier-Stokes applications. 7, 18, 54, 6 ' 

From a termwise point of view, a linearization of one- dimensional CDR equations can provide 
insight into the distinguishing characteristics of each term. Upon method-of-lines discretization using 
high-order, finite-difference techniques, the CDR equations may be written as a system of ordinary 
differential equations (ODEs) and analyzed with 

- X C U 4 \ D U • A'7 . (1) 

at 

where z c = A c (At), z D = X D (At), z R = \ R (A t), and (At) is the time step. The discretized convection 
term contributes scaled eigenvalues, z c , that are predominately imaginary while the diffusion terms 
have predominately real scaled eigenvalues, z D . Reaction rate eigenvalues are mostly real and may give 
rise to relatively large scaled eigenvalues, z R . Based on this knowledge, one might seek to construct a 
new method based on two separate methods: one optimized to smaller eigenvalues of the convection and 
diffusion terms and one that is capable of dealing with very large reaction-rate eigenvalues. It should 
be remarked that because of the high sound speeds introduced by the compressible equations, z G \ is 
generally larger than I* 0 !, e.g., in the DNS of a hypersonic boundary layer resolved down to y + = 1. 
Incompressible flows, governed by index-2 differential algebraic equations, generally have \z D \ > \z c \. 

If the stability domain of the integrator contains all values of z c , z D , and z R , then stable integration 
can be done. For accuracy purposes, the integration must proceed no faster than the fastest relevant 
physical processes contained within the governing equations. A situation may arise, however, where 
stable integration of the discretized governing equations can only proceed at a time scale substantially 
faster than any physically relevant time scale of the continuum-based compressible equations. This may 
render a numerical method unacceptably inefficient. It may occur in regions of intense grid clustering or 
while using stiff chemistry, but may be caused by other issues like interface boundary conditions within 
a multidomain formulation. There are two possible strategies to obtaining a solution at a reasonable 
cost when this occurs: change the governing equations and hence their characteristic time scales or 
change the numerical method. We choose the latter. 

Partitioned Runge-Kutta, methods may be designed to allow' for the partitioning of equations by 
term, gridpoint, or equation. Implicit-explicit (IMEX) partitioned methods developed to date for 
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first order ODEs have usually considered a partitioning based on terms. They have combined ex- 
plicit Runge-Kutta (ERK) schemes with variations on either the diagonally implicit Runge-Kutta 
(DIRK) 5 ’ 10 ’ 17 ’ 19 ’ 31 ’ 51 ’ 61 ’ 74 ’ 75 ’ 76 ’ 77 or Rosenbrock family of methods. 12 ’ 39 ’ 52 ’ 57 ’ 61 ’ 74 ’ 75 ’ 76 ’ 77 In 
this paper, methods are derived using stifhy-accurate, explicit, singly diagonally implicit Runge-Kutta 
(ESDIRK) schemes for their stability properties and higher stage-order. Of the methods currently 
available in the literature, some exhibit lower-order coupling errors, coupling stability problems, no 
error control, and poor ERK or DIRK/Rosenbrock stability properties. The new schemes endeavor to 
address all of these shortcomings without falling prey to new ones. 

A multistep 4 ’ 16, 35 approach to IMEX schemes is also possible. Higher-order SBDF methods given 
by Ascher et al., 4 based on BDF methods for the implicit portion, are of the form 


fc-i 


k - 1 / 


U(n+i) = _ J2 + (A t)fa Y, 


(-i yjb! 


7=0 


aU fc -i-0!0+l)! 


ri”piik + 


implicit ’ 


where a,- and fa are the values for the k-step, order-k, BDF methods. 45 These methods may be quite 
effective for equations whose stiff terms give rise to predominately real eigenvalues. At third-order and 
above the BDF methods are only A(a)-stable; they lack A-stability but still have a vanishing stability 
function for extremely stiff eigenvalues. Transient events like chemical ignition may be less appropriate 
for SBDF methods because they lack stepsize adjustment. 

The goal of this paper is to provide complete methods for solving CDR equations using additive 
Runge-Kutta methods. By this, it is meant that beyond having good accuracy and stability properties, 
there are high quality embedded methods, error controllers, dense output, stage-value predictors and 
implementation guidance. This is done in a two fold manner. First, a general overview is given 
of the coupling of N different Runge-Kutta methods for first-order ODEs whose right hand side is 
the summation of N terms: N-additive Runge-Kutta methods (ARKn). Second, matters are then 
specialized to the case of N = 2 : additive Runge-Kutta (ARK 2 ) methods using an implicit method 
that allows accurate and stable integration of the stiff terms while either integrating the nonstiff terms 
at the linear stability of the explicit method or integrating the entire method at some chosen error 
tolerance. 

We do not consider the topics of storage reduction, contractivity, dispersion and dissipation, reg- 
ularity, or boundary error. Sections 1 and 2 will provide an introduction to ARIvn methods and the 
motivation for their application to Navier- Stokes- type problems. Specialization of additive methods to 
ERK-ESDIRK combinations will be reviewed in section 3. Third , fourth-, and fifth-order schemes will 
be considered in sections 4, 5, and 6. Testing of the ARK 2 schemes on a one- dimensional, convection- 
diffusion-reaction test problem and three, two-equation, singular-perturbation problems is discussed in 
section 7. Merits of the additive schemes are discussed in section 8 and comparisons are made with 
existing Runge-Kutta and multistep IMEX methods. In section 9, conclusions are drawn as to utility 
of the various schemes. Appropriate appendices are also included. All new methods presented in this 
study were solved completely using Mathematica. ' 4 ’ ‘ 2 Coefficients of the new schemes are provided to 
at least 25 digits of accuracy. 


2. N-Additive Runge-Kutta Methods 

2.1. General 


Following Araujo et al., 3 ARKn methods are used to solve equations of the form 

N 

= F(U ) = £F [ " ] (E0, 

v = 1 


dU 

dt 


3 


( 2 ) 



where F(U) has been additively composed of N terms. They are applied as 

N s 

U (i) = u {n) + (At)J2J2 a[ if FW \ uU) ^ (3) 

i/=i j = 1 

Ns Ns 

u (n+ 1 ) = L T(n) + ({/«), #(»+ 1) = {tffi) + (A/) ({/«), (4) 

i/=l « = 1 i/=l 7 = 1 

where each of the N terms are integrated by its own s-stage Runge-Kutta method. Also, {/A) = U(t A) ), 
{/M = U (ti n ) + c;A{^ is the value of the {/-vector on the A 1 ' stage, and {/(” +1 ) = {/(/(") -f A{). Both 
f/M and {/( re+1 ) are of classical order q. The {/-vector associated with the embedded scheme, {/( ra+1 ), 
is of order p. Each of the respective Butcher coefficients oW, b^'\ and c[^, ;/ = 1,2,---, JV are 
constrained, at a minimum, by certain order of accuracy and stability considerations. 

2.1.1 Order Conditions 


Order-of-accuracy conditions for ARKn methods may be derived via N-trees. 3 These N-trees re- 
semble the traditional Butcher 1-trees 8, 26, 2 ' but each node may be any one of N varieties or colors. 
Expressions for the equations of condition associated with the c/* 1 - order N-trees are of the form 




a 


(?) _ 2!_ 

%] q\ 


1 

a 




where 


*K) — 




(5) 


where and 4> j are scalar sums of Butcher coefficient products and 1 < n < N q is used to 

distinguish between the many possible color variations of the k th 1-tree of order q. Both a and a are 
color dependent: i.e., for a given 1-tree, the many corresponding N-trees may have different values of 
a and cr depending on the details of the node colorings. Tree density, 7, is color independent and 
consequently so is the product of a and <j, a a = ql/'f. Order conditions t^ 3 \ t^, 7g 5 79 , r!f^ Q 16 17 lg 20 
given in appendix A, never exhibit color dependence. When an equation of condition, 77 j^- , is made to 

vanish, color dependence is immaterial because = I/7. For equations of condition that are not 
made to vanish, color dependence must be taken into consideration to accurately assess the leading 
order error terms. Order conditions for partitioned Runge-Kutta methods have also been derived by 
Jackiewicz and Vermiglio 3 ' following an approach of Albrecht. 


2.1.2. Coupling Conditions 


Aside from satisfying the order conditions specific to each of the elemental methods of the ARKn, 
one must also satisfy various coupling conditions. One may write the total number of order conditions 
for a general ARKn associated with each particular root node coloring, oq , using the expression 24, 26, 62 

OO OO AT [JV] 

n (1 - ■ m 

i= 1 i= 1 


At order i = {1,2, 3, 4, 5, 6, • ■ , a general ARKn method has a total of /Vo; order conditions, where 

<*! 1] = {1,1,2,4,9,20,---}, af ] = {1,2,7,26,107,458,---}, af ] = {1,3,15,82,495,3144,---}, a[ 4] = 
{1,4,26,188, 1499,12628, -■■}, and af ] = {1,5,40,360,3570,37476, ■}. Some of these, Naf\ are 

order conditions of the elemental methods which compose the ARKn- This implies that N(a^ — 0:[^ j 
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of the order conditions are composed of portions of different elemental methods. These are coupling 
conditions. As both order of accuracy and N increase, their numbers grow explosively. Table 1 shows 
their numbers for orders up to six and for each 1-tree listed in appendix A with which the coupling 
conditions are associated. Also contained in the table is the type of order condition: quadrature(Q), 
subquadrature(SQ ), extended subquadrature(ESQ), and nonlinear(NL). 68 It should be noted that there 
are no more than N * s(s + 1) independent Butcher coefficients to satisfy all of these order conditions. 

Obviously, without some type of simplifying assumptions, even fourth-order methods appear rather 
hopeless. The simplist remedy is to require all root or canopy nodes of the N-trees to be the same, or 
effectively colorless. In terms of Butcher coefficients, this is done by setting b]'^ = (root nodes) or 
c[^ = c\''^ (canopy nodes), = 1,2,---, N. Tables 2 and 3 show the results of these simplifications. 
Since there is always only one root node for any N-tree yet there may be many canopy nodes, assuming 
cf^ = c[A generally reduces the number of coupling conditions further than b]'^ = b \^ . Table 4 shows 
that matters become much more tractable when all root and canopy nodes are made equal. It may 
also be seen from this that as long as b\^ = b\^ and = e|^, an arbitrary number of independent 
third-order methods having the same number of stages may be coupled together with no associated 
coupling error. Tables 2 and 3 show that selecting 6-^ = b\^ or allows second-order error-free 

coupling of an arbitrary number of schemes. Choosing neither of these assumptions reduces error-free 
coupling, as seen from table 1, to first-order. 

2.1.3. Error 


Error in an elemental q th -order Runge-Kutta scheme contained within an ARKn method may be 
quantified in a general way by taking the principal error norm,' 11 


A ( g + i ) = || r ( ' i+1 >|| 2 



( 7 ) 


where rj ?+1 ^ are the error coefficients associated with order of accuracy q + 1. For embedded 
schemes where p = q — 1, additional definitions are useful such as 



fibp + 2 ) = 


^(p+2) 



z 


MS 


— A {p+ 1) = |if(p+ 1 ) 

P v 


(j(p+2) _ 


f(P+2) _ r (P+ 2 )|| 2 
i(P+l) 


£■(?+ 2) _ 


A(p+2) 

i(P+l)’ 


( 8 ) 

( 9 ) 


where the superscript circumflex denotes the values with respect to the embedded method. In the case 
of ARKn methods, we generalize the traditional expression for to 






N 


[1] 

a q+l n N,k,(q+ 1 ) 

E E ( 


k = 1 


1 


.(?+i)V 

fc[n] ) ’ 


(10) 


where we have used to denote the total number of ARKn order conditions arising from all 

variations of the l: 1,1 1-tree at order ( q + 1). For example, for general N-trees with N = 5, q = 5, and 
k = 15, table 1 gives n- 5 , 15.6 = 4370 + 5 = 4375. If root and/or canopy nodes have been set equal, then 
n N,k,( q +i) i s given in table 2, 3, or 4, whichever is appropriate. Since there is generally no reason to 
assume any particular order condition is more important than another, it is prudent to consider them 
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Table 1: General ARKn Coupling Conditions 
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NL 

18 

132 
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1G20 

N 3 (N 2 + l)/2! - N 

ESQ 

1G 
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1120 
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NL 

22 
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1870 
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14 

87 
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42 
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>) 
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11 
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30 

240 

1020 

3120 
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ESQ 

22 

159 

G3G 

1870 

N 4 (N + l)/2! - AT 

>) 
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NL 

18 

132 

540 

1G20 

N 3 (N 2 + 1 )/2! - N 

M 

Jtb 

u. 
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22 

22 

159 
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G3G 

G3G 

1870 
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N 4 (N + l)/2: - AT 
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SQ 

14 

87 

31G 

870 
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20 

ESQ 

30 

240 

1020 

3120 

N 5 - N 

ESQ 

ESQ 

30 

30 

240 

240 

1020 

1020 

3120 

3120 

N 5 - N 
N 5 - N 

SQ 

SQ 

22 

30 

159 

240 

G3G 

1020 

1870 

3120 

N 4 (N + l)/2: - N 
JV 5 - N 

q = 6 


418 

3084 

12548 

3737G 

N(N — 1)(129GAT 3 + 193GA r2 + 229GN + 237G)/5! 

q < 0 
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3G31 

14201 

41271 

N(N - 1)(129GA 3 + 25G1 A r2 + 3511A' + 4040)/5! 
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Table 3: ARKn Coupling Conditions with c[^ = c[^ 
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0 
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20 
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0 
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2 
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SQ 
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20 
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SQ 

0 

24 
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10 

3G 

84 
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12 

42 

90 
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0 

0 
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20 
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f 
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SQ 
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70 
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2 

G 

12 
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NL 

G 

24 

GO 
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2 

G 

12 

20 
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SQ 

2 

C 

12 

20 

N(N - 1) 

4 0) 
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ESQ 

G 

24 

GO 

120 

N 3 - N 

NL 

14 

78 

252 

G20 

N 4 - N 

>) 

ESQ 

G 

24 

GO 

120 

N 3 - N 

ESQ 

G 

24 

GO 

120 

N 3 - N 

NL 

10 

51 

15G 

370 

N 3 (N+ l)/2! - N 
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24 

GO 
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ESQ 

14 

78 

252 

G20 

N 4 - N 
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17, 

ESQ 

14 

78 

252 

G20 

N 4 - N 
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18, 

ESQ 

14 

78 

252 

G20 

N 4 - N 

>) 

19, 
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14 

78 

252 

G20 

N 4 - N 

>) 

20 

b Q 

30 

240 

1020 

3120 

N 6 - N 

q = o 


1G4 
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2940 

7580 

N(N - 1)(6N 3 + 39 N 2 + 87 N + 114)/3! 
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3540 

8870 

N(N — 1 )(GjV 3 + 45iV 2 -|- 120 N + 18G)/3! 
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Table 4: ARKn Coupling Conditions with and c[^ = c[^ 
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all. In certain special cases one may be able to rule out extended subquadrature, or nonlinear order 
conditions based on the structure of the ODE at hand. To evaluate the value of the error estimate, one 
may evaluate the quality parameters B 0+ 2 ), C^ p+2 \ and E^ p+2 ' 1 using each from each elemental 

method or using every 

All embedded schemes considered here are applied in local extrapolation mode. For a given order 
of accuracy, one strives to minimize Based on experience with q = p + 1 ERK pairs, B^ p+2 \ 

C( p+2 \ and EtP+V are ideally kept of order unity. Another term, 


D 


Max {Kf l\ b i IK \Mi |}> 


(li) 


is usually kept less than 20. 

Because of the large number of order conditions associated with the embedded scheme of an ARKn 
relative to any one of its elemental methods, we allow for the possibility that the order of the main and 
embedded methods differ by more than one, as is customary. 66 A priori quality criteria for q = p + 2 
pairs does not appear to have been derived in the context of first-order ODEs. 


2.1.4. Simplifying Assumptions 

Butcher 8, 27 row and column simplifying assumptions will be helpful in designing methods because 
they can reduce and simplify the order conditions. By comparing tables 1 through 4, one quickly 
surmises that higher-order methods effectively require the use of the assumptions = b\^ = b t and 
= c-^ = c;. Also, without identical root or canopy nodes, application of Butcher simplifying 
assumptions would become very awkward; therefore, simplifying assumptions are considered in the 
form 



M q — i i 

Y a ij C j = — 1 > = f/=l, 

U q 

(12) 


E^r 1 «!? = ^( !-«?). j = h-,s, q = 

i= 1 ^ 

(13) 


2.1.5. Stability 


The linear stability function for N-additive methods is considered using the equation 


N 


F(U) = Y, 


l/= 1 


from which it is determined that the stability function is 10 

P(zM,zW,---,zW) 






Det 

I~Eli( 

>1aM) 

i+Elil 

(«• 


Det 

I -Ell (zMAM). 



(14) 

(15) 

(16) 


where — b\ y \ I = 6 } j , z^ = A^Af. and e = {1, 1, - - - , 1}. Stability for the embedded 

method is considered using 


B(zM,z®.---,zM) = 


Det 

i -Eli ( 

zM AM) 

1 + Ell 1 



Det 

I -Eli (J? lAH) 



(17) 
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where bM = For nonlinear stability, which we do not pursue, one may consider 


N 


F(U) = £AH(F)F. 


( 18 ) 


i/=i 


Defining = iz^, z^\ ■ ■ ■ , 2 ^ j as the values of at times t) n ( + c,-Af, the Runge-Kutta K 

function is given by 


K{ = 

2.1.6. Conservation 


Det 


zMah] 

\ + Elr \ 

> [ 

"le® bH T | 

Det 

i-Eli {zWa(' 




(19) 


Conservation of certain integrals or invariants may also be of interest in additive Runge-Kutta 
methods. 3, 28 Similar to the algebraic stability matrix, one may define 




i % j 


3 3 l '3 


(20) 


ARKn methods conserve linear first integrals, in general, only if = 0, and conserve certain 

quadratic first integrals, in general, only if 6-^ — = 0 and = 0 , where i, j = 1, 2, • • • , s, v,l-i = 

1, 2, • ■ ■ , iV, v ^ p. Conservation of cubic invariants with Runge-Kutta methods is not possible. 9, 28 ’ 36 

2.1.7. Dense Output 

The purpose of dense output 26, 50 has traditionally been to allow high-order interpolation of the 
integration variables at any point, t( n > + #Af, inside the current integration step where 0 < 9 < 1. It 
may also be used, albeit more cautiously, for extrapolating integration variable values to enable better 
stage value guesses when one or more of the elemental methods is implicit. 49 ’ 55 For an ARKn method, 
it is accomplished as 


N 

I. 

V—\ 2 = 1 


U(t^ + 9At) = C/W + (Af)^^6* H (^)KM(f r W), b* [l/] (9) = b*\ l/] 9 j , = 1) = b [ - ] , (21) 


3= 1 


where p* is the lowest order of the interpolant on the interval 0 < 9 < 1. By construction, b*^\o = 
0) = 0. Order conditions, at order m, for the dense output method are given by 


I ^ 

*(m) _ - 

k[n ] g / -j * i,k[n\ 


a0 n 


ml 


( 22 ) 


Setting m = q and 9 = 1, we retrieve (5). As with the main and embedded formulae, one may write 
terms like A*( p * +1 ( = A< p * +1 \8) to access the accuracy of the dense output method. When used as 
an extrapolation device (8 >1), the stability function z™, ■ ■ - ,z^- N \9) must be considered, 6 


R*{z^\z( 2 \---,z^ n \9) = 


Det 

i-Elil 

[z( v 

■IaM) 

I+Elil 


3®b*M r (0)y 

Det 

I-Eli (»AH) 



(2.3) 


where b*M = b*^\ Throughout this paper, we will assume = b*^ = b*-. 
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3. Implicit-Explicit ARK 2 Methods 

Although one could consider a triple partitioning of the convection- diffusion-reaction equations by 
using two-explicit methods and one implicit method, the increase in coupling conditions from N = 2 
to N = 3 as well as the complication to the stability function warrants the former. Simplifying 
assumptions might facilitate solving extra order conditions, but at A r = 3 more trees of similar error 
result in higher principal error norms. Linear stability of the hybrid explicit method using a real-axis 
optimized method for diffusion and an imaginary- axis optimized method for convection will have many 
five- and six-node coupling tall trees with which to contend. The equations are therefore cast in the 
form 

^ = K ns + F s , (24) 

where F ns represents the non-stiff terms and F s represents the stiff terms, and we consider implicit- 
explicit ARK2 methods. ERK methods are used to integrate the non-stiff terms. Stiff terms are treated 
with ESDfRK methods. 2, 32, 42, 43 Coefficients for the ERK and ESD1RK methods will be distinguished 
by and o-^, respectively. ESDIRKs offer the advantages of allowing L-stability, stiff accuracy, and a 
stage-order of two. They differ from the more traditional SDIRK 1, 2 ‘ methods by having an explicit first 
stage. A consequence of allowing a stage-order of two is that algebraic stability becomes impossible. 64 
As we will always invoke b = /y, b\ E ^ = b-p = b t and = c\ ^ = c; in this paper, their 
superscripts are henceforth superfluous. In general, an ERK method has 5(5 + l)/2 degrees of freedom 
(DOF) available to satisfy all order and any other conditions, where 5 is the number of stages. An 
ESDfRK has (s 2 + s + 2) / 2 available DOF. Combining the two into an ARK2 scheme, if each b, and 
each c, are made equal, (2s — 1) DOF are lost, leaving (s 2 —5 + 2) DOF. A further assumption of 
a^J = bj, the stiffly- accurate assumption, reduces this to (s 2 — 25 + 2) and facilitates L-stability as 

well as forces the stiff part of f/( n +i) t.o be computed implicitly, ft is particularly useful in cases of 
singular-perturbation type problems and, when combined with I- stability, generally tolerates stiffness 
better than non-stifHy accurate L-stable methods. Incorporating stiff accuracy and a stage-order of 2 
into all of the ESDIRKs, 42 the IMEX ARK2 methods then take the form 


0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


0 

27 

27 

0 

0 

0 


0 

27 

7 

7 

0 

0 


0 

c 3 

[E] 

“31 

a [E] 

“32 

0 

0 


0 

C3 

“31 

a W 

“32 

7 

0 


0 

Cg — 1 

“5-1,1 

a m 

“5-1,2 

a [E] 

“5,2 

a [E] 

“5-1,3 

a [Ei 

“5,3 


0 

0 

C5-1 

a !/] 

“5-1,1 

a W 

“5-1,2 

a [I] 

“s — 1,3 


7 

0 

1 



[-E] 

“5,5-1 

0 

1 

61 

62 

&3 


b s - 1 

7 


bi 

b 2 



b s -i 

7 


6l 

b 2 

bs 


b s - 1 

7 


6l 

b 2 

h 


65-1 

b . 


61 

b 2 

b 3 


b s - 1 

bs, 


where 7 = (+, i = 2,3, ■ • - ,5 and should not be confused with the density of an N-tree, 7, discussed in 
section 2.1.1. 

To identify the schemes derived in this paper, a nomenclature similar to that originally devised 
by Dormand and Prince is followed. 41 Schemes will be named ARK q(p)sS[q so ]X , where q is the order 
of the main method, p is the order of the embedded method, s is the number of stages, S is some 
stability characterization of the method, q so is the stage order of the implicit method, and X is used 
for any other important characteristic of the method. For .5', we use L to denote an L-stable ARK2. 
Distinguishing between L-stable methods that are or are not stiffly accurate is important: hence, we 
use A' as S' A to denote stiffly accurate. 

3.1 Design 
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3.1.1 Accuracy 

Both ERK and ESDIR.K methods are subject to the -[1, 1, 2,4,9,20} order conditions for orders 
{1, 2, 3, 4, 5. 6). These order conditions are listed up to sixth-order for 1-trees and fifth-order for gen- 
eral 2-trees in appendix A. With the assumptions b\^ = = bi and c[^ = = c,-, there are 

{0,0,0,2,13,63} coupling conditions at the same orders (see section II.15 26 ). At fourth- and fifth- 
order, these coupling order conditions are shown as bicolored trees in figure 1. 




Further reduction of the number and complexity of order conditions is possible by using Butcher 
simplifying assumptions. Unfortunately, they may conflict with one another. For instance, applying 
assumptions and D^\l,j) gives rise to two inconsistent equations at j = s, 

6*7 = b s (l-c 3 ), 0 = b s (i-c s ). (25) 

This implies that either b s = 0, c s = 1, or c 8 = 1 — 7. The first is unacceptable while the remaining 
two are contradictory: hence, the column simplifying assumption can only be applied to one of the 
methods. In conjunction with the stiffly accurate assumption that forces c s = 1, it may only be used 
on the ERK. We avoid the option of setting both and all) to zero because of the complication that 
will arise in enforcing L-stability and the possibility of explicitly computed stage values. Additional 
interscheme conflict occurs upon imposition of £7^(3,*) and £7^(3, i). A stage order of two on the 
stiffly accurate ESDIRK is imposed by enforcing £7^(2,*) = 0 for i = 2, ...,(.s— 1). A stage order of 
three is impossible because of the second stage where a y c j — c !/3 = 4y 3 / 3 ^ 0. Reducing the 

truncation error of the second stage is clearly facilitated by smaller values of 7. 


3.1.2 Stability 


Linear stability for an ERK-ESDIRK ARK2 method with b\^ = b^ 
stability function 10 


b is analyzed using the 


R(zl E \zW) 


Det [i - zW AW - AW + (zW + e ® b T 
Det [I - iflAlfl] 


P(z^ E \z^) 


(26) 
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where 


P(z [E \z [I] ) 

Q(z [I] ) 


eIe^)']^)'. 

i = o [j=o ) 

i + E«(- J/1 )‘ = (1-7J' 1 ) 

i=l 


(27) 


^{((s-D-ky.kiJ 



(28) 


To recover the ERK method only, pij = poj and q t = 0. The ESDIRK is retrieved with p % j = p{Q. 
Both the ARK 2 and ESDIRK methods share the same q t coefficients. I 11 all cases poo = 1. A total of 
(5 — l)(s — 2)/2 of the pij coefficients are coupling stability terms. For an ESDIRK to be L-stable, it 
is required that 7 > 0 so that the stability function remains analytic in the complex left-half-plane, 
the method must be I-stable, and p s ,o = p s -i,o = 0 so that the stability function vanishes as z^ tends 
toward infinity. I-stability of the ESDIRK method is determined using the E- polynomial 2 ' given by 


E(y) 


Q(+iy)Q(-iy)- P{+iy)P(-iy ) 


s 


E Ppy'- ■ 

j=o 


(29) 


where i = \/—l and P(Eiy) is composed of only the p t 0 terms. Imaginary axis (Instability requires 
that E(y) > 0 for all real values of y. It is sufficient but not necessary to have all E 2 j > 0. An L-stable 
ESDIRK will have f? 2 s = 0. An order q ESDIRK will have E^j = 0 for 2 j < q. 

Above and beyond L-stability of the ESDIRK method, it may be useful to control the damping 
of the large scaled eigenvalues, zB\ at intermediate stages. 69 The internal stability function at the 
n th -stage may be constructed for DIRIv-type methods by using portions of the matrix. Denoting 
aW, i, j = 1 , 2 , • ■ -, n as A and , j = 1 , 2 , • • ■ , n as Bj, the internal stability function is given by 


Rt ] < 


c mt 


(* [/1 ) = 


Det \l - z^A [I] + z^S 0 B 1 


p( n ) ( , 

1 int 


f/]) 


Det [1 - zV\AW\ 


q 


w) 


(30) 


where X is the (n X n) identity matrix, and £ is the one- vector of length n. Our concern will be the 
value of ( — 00 ). P-^ is, in general, a polynomial of degree n— 1 in z^ because A n j = Bj while 
is, in general, of degree n. Consequently, SDIRKs have 00 ) = 0. ESDIRKs, with reduced 

to degree n — 1 because a = 0, do not generally satisfy — oo) = 0. 

In terms of step-wise stability, choosing the stiffly accurate assumption forces p s . 0 = 0. Placing 
cijj = 0 forces p s -i,o = 0 but sacrifices (5 — 1) DOF and the possibility of higher stage-order. A 

consequence of setting a ^ = 0, what effectively distinguishes the ESDIRK from the SDIRK, is that it 
forces q s = 0. Achieving an L- acceptable stability function for the ARK 2 , 


R{ 




:W) = 


1 + h {pa-l,l^ + Pa-1,0 } (^) + p 3 , 0 (^) 


1 + ■ ■ ■ + (- 7)' -1 (* W ) 


5 — 1 


(31) 


and not just the ESDIRK, is now more complicated because of u. I 11 both the ESDIRK and 
IMEX ARK 2 cases, p a , 0 and p s -i.o must vanish for L-stability, but the IMEX scheme must also satisfy 
Ps-1,1 = 0. Several of the methods given by Ascher et al. and Griepentrog 19 do not account for this 
and consequently have R{z'P^,— 00 ) depending linearly on z<- E 1. Similar comments apply to pij and 
P^j = X^ = i Pijk^i the coefficients of the P polynomials for the embedded and dense-output formulae. 
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A benefit of a zero-column in , i = 1 for an additive method is that on problems such as Kreiss’s 
problem, 13, 2 ' which may act like an index- 2 differential algebraic equation, the initial value of the 
algebraic variable is not propagated along with the solution. 38 Asclier et al. 5 and Calvo et al. 10 both 
choose a^j = bj as well as zero-padded SDIR.Ks = 0) and consequently their methods perform 
relatively well on Kreiss’s problem. 

3.1.3 Conservation 


As the scheme weights, bi, for the current ER.K and ESDIRK methods are the same, linear first 
integrals will be conserved. Certain quadratic first integrals will, however, only be conserved if linear 
first integrals are conserved and 

M t I] = b i a [ $ + b j af-b i b j = 0, (32) 

vanishes. We will adopt the point of view that although none of the methods will make this vanish, 
minimizing the magnitude of this matrix, 

n<u = (33) 

should enhance conservation characteristics. 


3.2 Implementation 


3.2.1 Stage Values 

Using the definitions F^) = F ns (u^\ t( n ) + c^-At) and Fg J ' 1 = F s (jJ^, + CjAt'j , one must 

solve 

C7(0 = L j{n) + x(0 + (A f) 7 .F s (‘>, i > 2, V') = (A t)J2 («[f F<£ + a^F^) , (34) 

3 = 1 


where X (0 is explicitly computed from existing data. Combining this with an appropriate starting 
guess, a modified 23 ’ 33 ’ 58 > 60 Newton iteration provides [00 and F$ \ In cases where direct methods 
are appropriate, this is accomplished by solving 



[/(»>) + + (A tfrF, (u^) , (35) 


where the subscript k denotes the value on the k th iteration, M = [I —j(At)(dF s /dU)\k\ is the iteration 
matrix, and d(0 = ([ 7(0 _ [/OA is the displacement. On the & th iteration one has 


Mdf = rf, + (36) 

where r*(0 is the residual. The iteration is terminated when either rijO (displacement test) or r[0 
(residual test) are sufficiently small, 33 ’ 65 

(*) (0 
^residual — 1 ” ( F ''' ' / ■ - ’ . ^displacement — 11 0 , C -S 0.005, 


(i) 

where e is the integration error tolerance and c is the tolerance ratio. F^s may then be computed using 

[7(0. Inexact 14 ’ 15, 30, 40 ’ 4 '’ 53 Newton methods may be more appropriate for larger systems of coupled 

equations. 


15 



3.2.2 Stage- Value Predictors 


Stage-value iteration convergence rates may be substantially improved and convergence failures 
may be minimized by choosing a good starting guess. The most primitive approach to obtaining a 
guess for the integration variables at the next stage is to use the most recent stage values; the trivial 
guess. An oftentimes better way to obtain stage-value starting guesses is by using a dense output 
formula. 27, 32, 49, 55 Second and later steps may use the function evaluations from the previous step to 
extrapolate into the current step. Stage- value guesses for the stage of the step n + 1 are obtained 
using function evaluations from step n as 


U®(t n + 9iAt) = U (n) 


0i = 1 + rci, 


+ (A + 

3=1 

(A t)( n+1 > 
r ~ (A f)(") 


(37) 

(38) 


Shortcomings of this approach include order- of- accuracy reduction when an interpolation formula is 
used in extrapolation mode and the introduction of instabilities into the extrapolated guess. As is 
sometimes done with the implicit error control estimate when substantial stiffness is present, one 
may wish to smooth the predicted stage value by multiplying it by the iteration matrix. 32 ’ 59 More 
sophisticated predictors have been derived for DIRK methods. 29, 44, 56 We do not consider these, in 
part, because computer memory management may become a problem. To conserve memory usage 
during extrapolation, all s-stages may be estimated at the beginning of the step and function values 
might then be overwritten by stage- value guesses, one equation at a time. For large r, the trivial guess 
may be more prudent. 


3.2.3. Error and Step-Size Control 


Step-size control is a means by which accuracy, iteration, and to a lesser extent stability are con- 
trolled. The choice of the (At) may be chosen from many criteria, among those are the (At) from the 
accuracy based step controller, the (At )i nv i sc ia and ( At) v ; gcou3 associated with the inviscid and viscous 
stability limits of the ERK, and the (At); t , er associated with iteration convergence. 23 If error control 
reliability is sufficient, CFL numbers may be removed as their function would be superfluous. For 
q = p + 1 pairs, one could consider timestep control of the IMEX schemes using I-, PI-, PID-, or 
PC- controllers 20 ’ 21 > 22 > 63 


(A t)[ n+1] 

(A4™ +1) 
(A ttm ] 

(At)fe +1) 



II^IU 

(A 

(At)*" -1 ) 


(39) 

(40) 

(41) 

(42) 


with k ss 0.9, t is a user specified tolerance, and p is the order-of-accuracy of the embedded method. 
The I-, PI, and PID-controllers are appropriate to explicit methods. Implicit methods use either I- or 
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PC- controllers. A PID- controller is considered because many stability optimized ERIv methods are SC- 
unstable with a Pi-controller for eigenvalues on the real axis. 41 Its characteristic roots are those of the 
equation () 3 + (pa — 1)<( 2 — p/3( + P7 = 0. Individual controller gains are obtained via kj = p(a — f3 + 7), 
kp = j)(/3 — 2q), and kp = p') . One may wish to keep controller gains fixed, independent of stepsize 
changes. This may be done by preselecting gains, then using u> n = and 


pa 


kj + kp + f ’ n } k D , 
\ J- + / 


p/3 = [k P + 2 uJ n k D ] , 


PI = 


\ u) r 


k D . 


(43) 


In the present context, we have selected kj = 0.25, kp = 0.14, kp = 0.10, or a = 0.49/p, f3 = 0.34/p, 
7 = 0.10/p when u> n = 1, giving characteristic roots of {—0.518,0.247,0.781}. The controller is SC- 
stable at all stability boundary points of the ERK in the complex left -half- plane for each of the three 
proposed methods. Inclusion of second- derivative gain, a PIDD 2 -controller, was not found to enhance 
control. The term is given as either 

^(re+l) _ jj(n+ 1) _ jj(n+ 1)^ fi(n+ 1) _ jy- 1 1) _ //(n+1)^ ^ (44) 

depending on whether sufficient stiffness is present to require smoothing. 32, 59 Both approaches may not 
behave well if substantial order-reduction occurs due to extreme stiffness because, at a minimum, p no 
longer reflects the actual order of the embedded method. Stability based time step limits involving the 
inviscid and viscous CFL numbers are given by ( A/ )i nv i sc id ~ A(Aa:)/(w+a) and (Af) v j scous \ v (Ax) 2 /v 
where a is the local speed of sound, u is the magnitude of local fluid convection speed, and v is an 
appropriate diffusivity of either mass, momentum, or energy. 41 For implicit-explicit methods, we select 

(Ai) ( " +1) = Min {( A/)^+J>(Af)i nviwid ,(A/) viscou ,,(A/) iterat ion} • (45) 


It remains to be determined which controller! s) are best suited to IMEX methods. In cases where 
q = p+ 2, we follow Tsitouras and Papakostas 66 and consider a modified I-controller. but not the PI-, 
PID-, and PC-controllers. In this case 


(Af)j” +1) = k(A/)<”> 


{/ 2 (Af)||^+ 1 || 30 


1 

P+1 

? 


where / 2 ~ 10 is experimentally optimized. 


(46) 


4. Third-Order Methods 

A third-order, 3(2) pair, ARK 2 scheme is designed in four stages using = b^\ = cj^, and 

simplifying assumption (7^(2, i). Stiff accuracy and a stage-order of two are incorporated into the 
four-stage ESDIRK method. The main method is obtained by solving for the .s 2 — 2s + 2 = 10 DOF 
(s = 4) using 

T [ k) = 0, k = 1,2,3, p 30 = P3 1 =0, c 3 = 3/5, , 

[i N 3) = 0, [£] ^i 4) = 1/35, E -=i <$ Cj = c\! 2 , % = 2, 3, 

with c 3 and being used for optimization. Total fourth-order principal error is AW = 0.07217. 

Linear stability of the ERK is given by (A,A„) = (1.24,0.92). Solving p 3 0 for 7 results in a cubic 
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equation, 67' 1 — I87 2 + 97 — 1 = 0, having three roots. Only one of these gives I-stability, 7 « 
0.435866521508458999416019. An embedded method is found by enforcing 

= 0, k = 1,2, p 4O = 0, i> 30 = — 3g 3 /40. (48) 

Second-order dense output is achieved by solving for b * = b*j0^ using 

2 

rl (k) = 0, k = 1, 2, = 64, = P.; 12 = 0. (49) 

i=i 

Its properties are summarized in table 5. Coefficients of the scheme, ARK3(2)4L[2]SA, are given in 


Table 5: ARK3(2)4L[2]SA Dense Output Method 


Property 

8 = 1 

8 = 2 

9 = 3 

8 = 4 

9 = 5 

A*( 3 )(6>) 

0 

1.098 

4.809 

12.87 

27.02 

A*W(0) 

0.07217 

2.280 

14.18 

48.23 

121.9 

R*(zW,~ oo,6) 

0 

5.789 

18.37 

37.74 

63.89 


appendix D. Characteristics of the method are listed in appendices B and C. 


5. Fourth-Order Methods 


A five-stage ARK2 method using a stiffly -accurate, L-stable, stage-order 2, ESDIRK method must 
satisfy 247 4 — 967 3 + 727 2 — I67+ 1 = 0. Of the four roots to this equation, only one leads to an L-stable 
method resulting in C2 = 27 % 1.14563212496426971081600277. Further, the minimum principal error 
associated with the two free-parameter family of five-stage, fourth-order, stage-order 2, stiffly-accurate, 
L-stable ESDIRKs is A* 5 ) = 0.03855, approximately 15 times greater than SDIRK4. 2, In spite of these 
shortcomings, one may construct a fourth-order, 4(2), ARK2 pair using identical root and canopy nodes 
as well as row' simplifying assumption C^(2, i). With 17 main and 5 embedded available DOF, one 
may solve 


rf * = 0, k = 1, 2, 3, 4, £j=i a^cj 
P40 = P4r = = [ ^ 4) = 0, 

c 3 = 50/100, C4 = 95/100, ff k) = 0, k = 1,2, 


- C 1 l^i * — 2, 3, 4, 

1 c k = 1/4!, 

b^ 3) = #50 = 0 . #40 = -?4/10, 


(50) 


and obtain a leading order principal error norm of A( s ) = 0.07664. Linear stability limits for the explicit 
method, in terms of the inviscid and viscous CFL numbers, are (A,A„) = (1.38,0.67). 

For the design of a fourth-order, 4(3), ARK2 pair, we again use the simplifying assumptions Vp = 
bp and c-^ = c-^. Only stiffly accurate, stage-order 2, ESDIRK methods are employed. Using six 
stages permits s 2 — 2s + 2 = 26 DOF (s = 6). The value of 7 must be 


0.2479946362127474551679910 < 7 < 0.6760423932262813288723863, (51) 

for I-stability. 2, Besides facilitating better iterative convergence of the modified Newton method, smaller 
values of 7 tend to result in low r er truncation error in the implicit method. For simplicity, we use 
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7 = 1/4. With assumptions (7^(2, i) and (7^(2. i), ARK4(3)6L[2]SA satisfies 


r W _ 


= 0, k = 1,2, 3, 4, E U a ih = c */ 2 > * = 2 ’ 3 ’ 4 ’ 5 > Ei=! «lj C 1 = Cf /2, * = 3, 4, 5, 6 

1. — n _ _ [E]_( 4 ) _ [J] _( 4 ) _ IB _ /,.„[£] _ n 

°2 — P.50 — F51 — r 3 — r 3 — Z_i = l °»®j2 — l^i= 1 C a i2 ~ 0, 

c 3 = 332/1000, c 4 = 62/100, c 5 = 85/100, 7 = 1/4, 

[ e V 5 (5) = 1/8000, = 1/135, = 1/1250, 


( 52 ) 


where = 0.01224 and (A,A„) = (2.01,1.06). The embedded method for this scheme is found by 
solving 


T] h) = 0, k = 1,2,3, b- 2 = p 6 0 = 0, p 50 = -3g.5/20. 


(53) 


Dense output may be approached by either maximizing accuracy or stability. If the method is to be 
used for interpolation, then a third-order method is appropriate. For extrapolation, stability is more 
important and a second-order method is constructed. Third-order dense, or continuous, output is 
achieved by solving for b* = Ej=i b* ;j 9 :i with 

r* (k) = 0, k = 1, 2, 3, b * = p* 0 = 0, p* 13 = -3 95 , pi 02 = -3 q 5 , E?=i Kj = h- (54) 

Second-order dense output is achieved with p* = 2 by solving 

T P ^ = 0; k = 1,2, &2 = PsO = Ps 1 = P502 = 0, Ej=l ^6 j = ^6- (55) 

Properties of these two methods are given in table 6. Characteristics and coefficients of the scheme 


Table 6: ARK4(3)6L[2]SA Dense Output Method 


Property 

6 = 1 

II 

to 

9 = 3 

II 

4^ 

II 

2 nd - Order 






A*^(0) 

0 

1.106 

5.049 

13.56 

28.38 

A*( 4 \6) 

0 

2.459 

14.79 

49.52 

124.1 

R*(z^,-oo,9) 

0 

-1 

-2 

-3 

-4 

3 rd - Order 






AA 4 )(0) 

0 

2.206 

12.86 

42.48 

106.5 

A*( 5 \0) 

0.01224 

3.852 

30.61 

132.1 

410.2 

oo,0) 

0 

21.2 - 10.4 zW 

92.6 - 49.1 z^ 

242 - 134 

499 - 284 z^ 


ARK4(3)6L[2]SA are given in appendices B, C, and D. 

6. Fifth-Order Methods 

In principle, one may construct a fifth-order ARK2 method in seven stages using simplifying as- 
sumptions $p = bp\ cp = cp\ (7^(2, i), D^(l.j’), and (7^(2, i). L-stable methods require that 

0.1839146536751751632321436 < 7 < 0.3341423670680504359540301. (56) 
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With the 37 DOF available, minimally, the following 33 equations must be solved 


r (*) _ 


= 0, k = 1, 2, • • ■ , 5, E-=i b t af = 6,(1 - Cj ), j = 2, 3, ■ • • , 6, 
ELi “ [ ijCj = cf/ 2 , i = 2, 3, • • • , 6, £y=i a [ p ] Cj = c?/2, i = 3, 4, • • • , 6, 


/, _ „ _ „ _ [_B]_( 5 ) _ I7]_( 4 ) _ [Z] _ ( 5 ) _ 1 J „2 I /rj _ n 

02 - Peo - Pei - 1 - 1 J T 3 - l J r 4 5]8 _ ^ E,-,j=i Eu- E “ V 51 - u , 




ELi = E;L, b iCi a\ f 1 = E-=1 ft.-c.-aS5 = Ei,j=i = Em=! b t a [ Vaf 2 = 0. 


(57) 


To solve this system of equations, one must solve for at least one abscissa directly. Given the size of 
the system, it is more fruitful to temporarily ignore E;=i ft." c .-aS5 = 0 and E; j=i b i a p a p = 0? include 

^ r 20^ = specify C3, cq, C5, and Cq, and investigate if the scheme merits further effort. With these 
changes, both implicit and explicit methods are fifth-order but the coupling method is only fourth-order. 
We have been unable to find any promising solutions to this method. 

Adding a stage, we now consider an eight-stage 5(4) pair. Eight stages permit s 2 — 2s + 2 = 50 DOF 
in the main method and 8 in the embedded method. The primary difficulty in designing a 5(4) pair is 
reducing the number of embedded order conditions while simultaneously keeping the implicit portion In- 
stable and keeping the tall trees of the explicit method well placed. We select simplifying assumptions 
bP = bp\ cP = cf\ C^(3,«), D^ e \ l.j), and C^ E \'2,i). I11 addition, we set the lower elements 
within the second column of each matrix to zero. ARK5(4)8L[2]SA is constructed according to the 
following conditions 


rf } = 0, * = 1,2, •■■,5, E‘ =1 fti«!? = ft,(l 


ELiO-pCj = cf/ 2, i = 2, 3, - - - , 7, 


Cj), j = 3,4, ■••,7, 


Ej=i ajf cj = cf/ 2 , i = 3,4 r 


-,7, 


Ej=i up 0 ] = cf/3, i = 3, 4, • • • , 7, 4? = 0, i = 4, 5, • • ■ , 7, a[f = 0, i = 4, 5, • ■ • , 

h = h = p 70 = pn = [, ' ; M 5) = [/] T5 5) = E?=i b i a P = 5 Eij=i 4? c * - !/5! = 0, 


(58) 


' r ( k ) _ 


0, * = 1, 2, • ■ ■ , 4, b- 


= h = Pso = pj 4) = 0, Pro/g? = 1/5, 


c 5 = 92/100, cq = 24/100, c 7 = 60/100, 7 = 41/200 , 4? = -1/8, 


J£] 

l 7 6 


= -1/8, 


where A/ 6 ) = 0.006988 but (A,A„) = (0.43,0.96). Better 5(4) pairs probably require nine-stages. 

A fourth-order dense output formula is not possible. Third-order continuous output, however, is 
achieved by solving for b* = Ej=i b *j^ with 


E?=i b *c 


k - 1 


1/*!, k = 1,2,3, 


^3 — Pso — Pri — P703 — P702 — 0, Ef=i b 8j — (59) 


For ARK5(4)8L[2]SA, the dense output method is characterized in table 7. Coefficients of the scheme 


Table 7: ARK5(4)8L[2]SA Dense Output Method 


Property 

9 = 1 

9 = 2 

II 

CO 

II 

rO 

II 

A*W{0) 

0 

0.295 

4.434 

21.01 

63.59 

A*( 5 )(0) 

0 

1.581 

21.03 

107.8 

361.2 

R*{z^ e \ —00, 6) 

0 

-1 

-2 

3 

-4 


AR.K5(4)8L[2]SA are given in appendix D while method properties are given in appendices B and C. 

7. Test Problems 
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To test the IMEX ARK 2 methods that have just been presented, four separate test problems having 
adjustable stiffness are considered. Ultimately, all equations are singular perturbation problems 25, ' 3 
and each may be evaluated with either unperturbed or the more troublesome perturbed initial conditions. 2 ' 

7.1. Kap’s Problem 

Dekker and Verwer 13 investigate a nonlinear problem (experiment 7.5.2) originally given by Kaps, 

yi(t) = - (s' 1 + 2 ) yi (t) + £ y 2 (t ) = yi(t) - y 2 {t) - y%(t), (60) 

where 0 < t < 1 and whose exact solution is y\(t) = t/|(f), y 2 (t) = exp (— /). Equilibrium (unper- 
turbed) initial conditions are given by 2/1 ( 0 ) = 2 / 2 ( 0 ) = 1. The equations exhibit increasing stiffness as 
6 — » 0 and, in the limit of e = 0, the system becomes an index-1 differential algebraic equation system. 
This may be easily seen by multiplying the first equation by £ to obtain £?/i(t) = — (1 — 2 s) yi(t)-\-y 2 (t)- 
Upon setting £ = 0, it reduces to the algebraic equation yi(t) = y 2 {t)- I 11 an IMEX formulation, terms 
multiplied by £ _1 are integrated implicitly while all other terms are integrated explicitly. 

7.2. Van der Pol’s Equation 

Van der Pol’s (vdP) equation is an equation describing nonlinear oscillations where the solutions are 
damped (amplified) for large (small) values of iji , 26 ' 2 ‘ 

yi(t) = y 2 {t), y 2 (t) = £ _1 ((1 - yi{tf)y 2 (t) - t/i(t)) . (61) 

Unperturbed initial conditions are given by r/i ( 0 ) = 2, y 2 (0) = —0.6666654321121172. For partitioned 

integration, the first equation is integrated explicitly while the second is integrated implicitly. Van der 
Pol’s equation develops a very challenging boundary layer at time T ss 0.8 based on these initial 
conditions. Two test cases are chosen involving different time intervals: 1) 0 < t < 0.5 and 2) 

0 < t < 1.5. The first is used to study order reduction while the second tests the error prediction 
capabilities of the schemes and of the robustness of the error controllers. 

7.3. Pareschi and Russo’s Problem 

Pareshi and Russo 51 have constructed a simple test equation which contains both stiff and nonstiff 
terms, 


m(t) = -V 2 (t), = yi(t ) + e 1 (sin( 2 /i(<)) - y 2 (t)) . (62) 

Partitioning for an IMEX scheme, terms multiplied by £ -1 are integrated with the implicit method 
while other terms are integrated explicitly. Initial conditions may r be considered in two different forms. 
Equilibrium initial conditions remove any contribution of the stiff term in the initial conditions. This 
is accomplished with ?/i(0) = 7t/2, 2 / 2 ( 0 ) = 1. Nonequilibrium, or perturbed data is specified by 

replacing the condition on y 2 with 2 / 2 ( 0 ) = 1/2. 

7.4. One-Dimensional Convection-Diffusion-Reaction Problem 

A simplified one- dimensional version of the gas-phase, multicomponent, compressible Navier-Stokes 
equations with chemical reaction is tested. The simplification assumes no bulk viscosity, no thermal 
diffusion or its cross effect, 110 spatial gradients of the transport coefficients (p., A, pH/), no barodiffusion, 
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ordinary diffusion is representable by an effective Fickian diffusion coefficient, and no body forces are 
present. With these assumptions, one must solve the system 



P 


pu 


0 


" 0 " 

d 

pu 

d 

pu 2 -|- p 

+ 

4/l 1 d 2 u 
3 dx 2 

+ 

0 

dt 

pco 

dx 

(pc 0 +p)u 

*0 + ^0 + fit ) 2 + tirr=\pD,h,^) 

0 


_ pY t _ 


puY, 


L pD^ J 


Cji 


(63) 


where the three righthand side terms are Fq, F b, and Fr with the subscripts C, D, and R denoting 
convection, diffusion, and reaction, respectively. Also, u is the fluid velocity, T is the temperature, 
Yi are the species mass fraction, i is the species index that runs from 1 to lies (number of chemical 
species), p is the fluid density, p is the pressure, t is time, x is the spatial direction, eo is the total 
specific internal energy, p is the molecular viscosity, A is the thermal conductivity, I) : is the effective 
Fickian diffusion coefficient, hi is the partial specific enthalpy of species i, and <!;,■ is the reaction rate 
of species i. We consider the reaction rate only in modified Arrhenius form without pressure correction 
terms. Supplementary relations that are needed to solve this system are given by 


n j 11 C.O 

e 0 = u 2 /2 - p/p + ^ h,Yi, hi =h/+ c pi dT , e 8 = h t - R°T/W t , c p = ^ c pi Yi 

8 = 1 •' T ief i = l 

ncs / ncs \ — 1 

£ X = !» P = pRT = P R°T/W , W = Y i/ W i > ?p-c v = R = R°/W, 

8 = 1 \i = l / 

% = A k T th exp (Jpf) ^ K k = k f Jk eqk , <j> h = ^ , 

ncr 

*.■ = ^ e ? d - <r) 


k=i 


% n i § 


react 


prod -| 




(64) 

(65) 

( 66 ) 

(67) 


where fl° d are the stoichiometric coefficients for the products, are the stoichiometric coefficients 

for the reactants, rp^ are the third-body efficiencies for the reactants, a k is the third-body efficiency 
exponent, W is the average molecular weight, W t is the average molecular weight of species i. A k is the 
temperature prefactor for reaction k, j3 k is the temperature exponent for reaction k, E k is the activation 
energy for reaction k , R is the gas constant, R° is the universal gas constant, c p! - is the partial specific 
heat capacity at constant pressure, e; is the partial specific internal energy, kf k is the forward specific 
reaction rate constant of reaction k , k Tk is the reverse specific reaction rate constant of reaction k , k eqk 
is the equilibrium constant of reaction k, h9 is the reference partial specific enthalpy of species i , and 
ncr is the number of chemical reaction steps. 

This constitutes (ncs + 2)*(nx) equations that must be solved where nx is the number of grid points. 
In this work both convection and diffusion are deemed non-stiff, F ns = Fq, + Fr, and are integrated 
using the ERK method while reaction terms are treated as stiff, F s = Fr and are integrated using the 
ESDIRK method. To solve this system, a Jacobian of the stiff function with respect to the integration 
variables is needed. This may be done numerically or analytically. It is useful to recast the Jacobian 
as 


8F S _ dF s dV 
W ~ WdU' 

where F {p, pu, pc o, pi l, pi 2? ' ‘ ‘ 8 Pi (ncs— 1)1 8 I i/F ^ ‘ ^ ■ I 18 I 28 ‘ ‘ * 8 I (ncs— 1)} 8 

F s -(0, 0, 0, uq, U?2i ' ' * , d^ ncs _i) } • 


( 68 ) 
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For the purposes of this test problem, five species are included: H 2 , O 2 , OH, H 2 O, and N 2 . U T 
is then {p, pu, pe^, pY±{ 2 , pYo 2 , pYon 2 , pYn 2 o}. Wherever possible, an attempt is made to match the 
thermophysical properties of these molecules. A two-step reaction mechanism is employed having one 
reversible and one irreversible step 

H 2 + 0 2 r20H, ^H2 + 0H7 h 2 0, (69) 

, 2 k-i 

k 2 s 

where = ^ > k%. Values of specific reaction rate constants, fc,-, used in this test problem are not 
those of the supplementary relations above but have fixed AkT^ k and activation energy. In this way, one 
may simply adjust stiffness via the ratio & 2 /& 3 5 yet retain temperature dependence for the purposes 
of ignition. A parametric study identified the maximum stiffness (£ 2 /^ 3 ) for which the convection- 
diffusion-reaction system is stable with explicit time advancement. This value of fc 2 / U 3 is defined as a 
stiffness of e = 10°. A stiffness of 1(F implies that e _1 = fc 2 /&3 is HP times larger than its baseline 
stiffness value. Two levels of heat release were used in the study. The isoenthalpic case assumed that 
the enthalpy of formation of all the species was identical. The exothermic case assumed H 2 , 0 2 , N 2 
and OH to be identical, but H 2 0 was adjusted to yield approximately realistic flame temperatures for 
a hydrogen-air system. Derivatives are evaluated using sixth-order explicit stencils on a grid having 
401 points. Approximately 20 grid points are contained within one shock thickness. Spatial boundary 
conditions for the integration variables were specified using supersonic Euler conditions at the inflow 
and extrapolation conditions at the outflow. As gradients of flow variables were extremely small at the 
boundaries, specification error of conditions was deemed negligible. 

Initial condition for the hydrodynamic variables, {p, u,T}, are specified by using a precomputed 
normal-shock profile of air travelling at Mach 5. Each variable is nondimensionalized by its upstream 
value and integration is performed on the vector U = {p, pu, pea, pY{}. The nonequilibrium aspect of 
the initial condition consists of specifying a constant spatial distribution for species mass fractions. 
Stoichiometric reactant mass fractions are used. Initial values for reaction products are zero. Upon 
starting the simulation, isoenthalpic or exothermic chemical reactions are abruptly activated. Inte- 
grating the reacting shock wave through the spatial domain ten times with this species profile rapidly 
results in a consumption of 0 2 and H 2 behind the shock wave, an increase in H 2 0 behind the shock 
wave, and a small region of high OH concentration just behind the shock wave. This new profile for 
all integration variables is the equilibrium profile. Although testing integration methods with nonequi- 
librium initial conditions may seem somewhat contrived, fluid dynamicists rarely know their initial 
condition exactly. Oftentimes, the initial condition amounts to an educated guess, particularly inside 
a flame. 


8. Discussion 

Implicit -explicit, additive R.unge-Kutta methods from third to fifth-order are presented that allow 
for integration of stiff terms by an L-stable, stiffly-accurate, stage-order two, ESDIRK method while 
the nonstiff terms are integrated with a traditional ERK. Satisfied order conditions expressing splitting 
error are of equal order to those of the two elemental methods. Both the ESDIRK and ARK 2 methods 
have vanishing stability functions for very large values of the stiff scaled eigenvalue, — ► — 00 . Error 

control using a P ID -controller and dense output for interpolation and extrapolation are also provided 
for the new methods, unlike most existing methods. All constructed methods retain high stability 
efficiency in the absence of stiffness, — > 0. Each has been optimized to minimize the leading order 
ARK 2 error terms, minimize the size of the Butcher coefficients, maximize the stability envelope of the 
ERK, and maximize the conservation properties with respect to first integrals. The methods permit 
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partitioning of the ODE system by term, grid point, or equation. Numerical tests of the new schemes are 
conducted on a chemical reaction inducing propagating shock wave and three two-equation singularly- 
perturbed initial- value problems. The performance of these methods are compared to many existing 
ARK 2 methods: the (1,2,1), (2,2,2), (2,3,2). (2,3,3), (3,4,3), and (4,4,3) methods of Asclier et al., 5 
LIRK3 and LIRK4 due to Calvo et al., 10 a five-stage, 3(2) pair of Fritzen and Wittekindt (FW53), 1 ' 
ASIRK-3A from Shen and Zhong, 61 the LSSIRK-3A and LSSIRK-4A methods of Yoli and Zhong,' 5 and 
ASIRK-3A by Zhong 7 ' as well as the SBDF methods of Asclier et al. 4 Tests are conducted to determine 
stiffness leakage, efficiency, order reduction, error control quality, and dense output performance. No 
attempt is made to assess conservation properties. Characteristics of the various ARK 2 are listed in 
appendices B and C. 

8.1 Stiffness Leakage 

An essential requirement for the viability of stiff/nonstiff IMEX schemes is that the stiffness remains 
truely separable. If this were not the case then stiffness would leak out of the stiff terms and stiffen 
the nonstiff terms. It would manifest itself as a loss in stability or a forced reduction in stepsize of the 
nonstiff terms. A more expensive fully implicit approach might then be required, and hence, methods 
that leak substantial stiffness might best be avoided. We test for leakage on the reacting shock wave 
problem. There are two primary affronts that can be made to the integrator on this problem. The first 
is simply a very stiff reaction rate describing isoenthalpic or exothermic reactions. Secondly, one could 
provide an initial condition to the flowfield that is substantially different from the quasi steady-state 
solution. This nonequilibrium initial condition is accompanied by a strong equilibration process of the 
flowfield during the initial time steps. Ideally, this initial perturbation is damped during subsequent 
time steps. 

Thirty-seven existing and new IMEX ARK 2 arid three SBDF schemes are considered for testing. 
A mild test for stiffness leakage is to provide the integrator with an equilibrium initial condition and 
an exothermic reaction rate. Time steps are specified as CFL numbers where C’FL = u(Af)/(Aa:) 
and the velocity is Mach 5. Leakage generally falls into three catagories: insignificant, moderate, and 
catastrophic. Table 8 shows that for several methods the decrease in stepsize for the nonstiff method 
is in direct proportion to the increase in reaction rate stiffness. This constitutes, in our opinion, 


Method 

IflSW 



£ = 10 6 

Zhong, ASIRK-3A 

0.51 

0.01 

0.0001 


Yoh, SIRK-3A 

0.51 

0.01 

0.0001 


Shen, ASIRK-3A 

0.57 

0.05 

0.0005 


Yoh LSSIRK-3A 

0.52 

0.05 

0.0005 


Yoh LSSIRK-4A 

0.36 

0.06 

0.0006 



Table 8: Examples of catastrophic leakage. Maximum CFL as a function of reaction rate stiffness using 
equilibrium initial conditions. 


catastrophic leakage and a failure of the methods. They do not possess sufficient stability to be useful 
in the contexts that they might reasonably be expected to apply. We do not consider these methods 
further. 

A more severe test of leakage is a nonequilibrium initial condition with isoenthalpic reactions. In 
this case, table 9 shows that the two methods of Griepentrog 19 can be broken in this rather severe 
environment. Although much less severe than the leakage displayed in table 8, both methods of 
Griepentrog may be inappropriate for stiff computations. More reluctantly than above, we do not 
further consider these methods. It is interesting to ask why these methods have failed while other 
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Method 

O 

O 

rH 

II 

OJ 

£ = 10“ 2 

£ = 10“ 4 

to 

1 

o 

II 

OJ 

£ = 10“ 8 

Griepentrog (3-stage) 

0.37 

0.37 

0.37 

0.0007 

0.0 

Griepentrog (4-stage) 

0.44 

0.44 

0.05 

0.001 

0.0 


Table 9: Maximum CFL as a function of reaction rate stiffness using nonequilibrium initial conditions 
and a nonexothermic reaction rate. 


methods have not. One may inspect the internal stability function magnitude of the implicit method 
at stage i for — ■> — oo; R\'J t ( — oc). Unlike SDIRK methods for which /?j n j ( — oc) = 0, ESDIRK 

methods generally have nonzero values of oo). Table 10 shows that the final stages of both 

Griepentrog’s methods are noticably unstable. The IMEX Runge-Kutta methods of Ascher et al.. 5 


Stage 

1 

2 

3 

4 

Step 

Griepentrog (3-stage) 
Griepentrog (4-stage) 

+1.00 

+1.00 

+0.366 

+0.235 

-2.464 
+ 1.068 

-4.909 

-0.732 

+0.000 


Table 10: — oc) for the implicit methods of Griepentrog. 


Calvo et al., 10 Fritzen and Wittekindt, 1 ' and the present methods exhibited little to no leakage on 
either of these problems. If internal stability of the implicit method is a principal contributor to the 
breakdown of these methods, it is not surprising that the zero-padded A-matrices found in the methods 
of Ascher et al. 5 and Calvo et al. 10 do not leak stiffness badly because R ■„[( — oo) = 0 for all stages. 
The primary concession of this approach is that the implicit method cannot have a stage-order of two 
because <t 2 i i 1 a 22 - 

Our most severe leakage test for the reacting shock wave problem is to use a nonequilibrium initial 
condition and exothermic chemistry. This initial condition is severe enough to cause stiffness leakage 
for all of the Runge-Kutta methods when used in a fixed stepsize mode. Table 11 documents the 
progressive failure of several methods as stiffness is increased. Comparing this to the internal stability 


Method 

£ = 10° 

£ = 10" 1 

£ = 10“ 2 

£ = 10~ 3 

£ = 10~ 4 

£ = 10~ 5 

£ = 10~ 6 

£ = 10~ Y 

Ascher (2,3,3) 

0.38 

0.22 

0.17 

0.16 

0.12 

0.07 

0.01 

0.001 

ARK3(2)4L[2]SA 

0.57 

0.28 

0.20 

0.19 

0.14 

0.07 

0.06 

0.001 

Calvo LIRK4 

0.44 

0.44 

0.44 

0.44 

0.40 

0.15 

0.03 

0.008 

ARK4(3)6L[2]SA 

0.67 

0.49 

0.36 

0.34 

0.33 

0.16 

0.01 

0.001 

ARK5(4)8L[2]SA 

0.43 

0.41 

0.10 

0.05 

0.03 

0.001 

0.00 

0.000 


Table 11: Maximum CFL as a function of reaction rate stiffness using nonequilibrium initial conditions 
and exothermic chemistry. 


characteristics of the implicit methods in table 12, one may surmise that the ARK5(4)8L[2]SA is failing, 
in part, because of marginal damping at each stage. It is also interesting that ARK3('2)4L[2]SA and 
Ascher (2,3,3) wither similarly. One has strong damping at each stage while the other has strong 
damping at the end of the step. LIRK4, with strong damping at each stage and the step, arguably, 
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Stage 

1 

2 

3 

4 

5 

6 

7 

Step 

Ascher (2,3,3) 

+0.000 

+0.000 

- 

- 

- 

- 

- 

-0.732 

ARK3(2)4L[2]SA 

+1.000 

- 1.000 

-0.806 

- 

- 

- 

- 

+0.000 

Calvo LIRK4 

+0.000 

+0.000 

+0.000 

+0.000 

+0.000 

- 

- 

+0.000 

ARK4(3)6L[2]SA 

+1.000 

- 1.000 

-0.774 

-0.083 

-0.157 

- 

- 

+0.000 

ARK5(4)8L[2]SA 

+1.000 

- 1.000 

-0.732 

-0.649 

+0.856 

-0.967 

-0.353 

+0.000 


Table 12: ll[‘^ ( — 00 ) for IMEX implicit methods. 


appears to endure the best. The only IMEX methods that we are aware of that can integrate this 
problem in constant stepsize mode are the SBDF methods of Ascher et al. 5 In table 13, they show 
insignificant stiffness leakage on our most severe case. Notice that as the order of accuracy increases, 


Method 

s = 10° 

£ = 10~ 2 

£ = 10~ 4 

£ = 10~ 6 

SBDF2 

0.20 

0.17 

0.10 

0.09 

SBDF3 

0.14 

0.13 

0.10 

0.09 

SBDF4 

0.11 

0.10 

0.10 

0.09 


Table 13: Maximum CFL of SBDF methods as a function of reaction rate stiffness using nonequilibrium 
initial conditions and exothermic chemistry. 


the relative leakage decreases. If ARK3(2)4L[2]SA, ARK4(3)6L[2]SA, and ARI\5(4)8L[2]SA are used 
in conjunction with stepsize control, they are able to navigate through the strong initial transient of 
this problem. Although we cannot say conclusively why the SBDF methods remain stable while all of 
the IMEX AR.K 2 did not, it seems that both the order and stability matter. Runge-Kutta schemes 
may be able to satisfy R[^ t (—oo) = 0 at each stage but always have an overall stage order of one. The 
SBDF methods are L(oi)-stable to stiff eigenvalues but each step value is of design order. 

In practice, simulations of chemical systems often use exothermic reaction mechanisms and mod- 
erately nonequilibrium initial conditions. Once the computation is under way, each step would likely 
begin with nearly unperturbed initial conditions. For DNS of hydrocarbon flames using a compress- 
ible NSE formulation, spatial grid spacings may be of order ten microns at atmospheric pressure. 11 
Corresponding convective time step limits based on « 1 are of order ten nanoseconds. Under 
these conditions, detailed methane-air chemical mechanisms do not introduce time scales appreciably 
faster than those dictated by convective stability, implying e ~ 10°. Larger hydrocarbon mechanisms 
such as heptane-air may introduce timescales of order one femtosecond. Choosing to integrate at the 
convective limit, the stiffness would be approximately e k 10“ ' . 

8.2 Accuracy and Efficiency 

Beyond avoiding stiffness leakage, one would like accurate and efficient methods. Although our 
focus is principally on methods of third-order accuracy and higher, we will occasionally use a first-order 
method, Ascher (1,2, 1), and a second-order method, Ascher (2,3,2), for comparison purposes. These 
are chosen because their explicit methods have both non- vanishing CFL and viscous CFL numbers. At 
third-order, the first matter is to verify the accuracy of the methods: ARK3(2)4L[2]SA, Ascher (2, 3, 3), 
Ascher (3,4,3), Ascher (4,4,3), Calvo LIRK3, and Fritzen FW53. From appendix B, all methods are 
formally third-order. Ascher (4,4,3) and Fritzen FW53 use five-stages overall, four of them implicit, 
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while others use one or two less; hence, they might be expected to be relatively less efficient. Using 
only two implicit stages, Ascher (2,3,3) might be expected to be quite efficient. Ascher (2,3,3) and 
Ascher (3,4,3) have ARK 2 stability functions that depend on z^ E 1. All schemes of Calvo et al. and 
Ascher et al. have vanishing internal stability functions for the implicit method when > — 00 . 

Only ARK3(2)4L[2]SA and Fritzen FW53 have embedded methods, only Ascher (2,3,3) is not stiffly 
accurate, and only ARK3(2)4L[2]SA uses a stage-order two implicit method. Fritzen FW53 may not 
conserve linear first integrals well. Finally, values of 7 range from approximately 0.4359 to 1.0. 

Accuracy and efficiency tests are conducted using equilibrium initial conditions and exothermic 
chemistry on the CDR problem. All methods exhibit third-order accuracy in the absence of stiffness 
where error is given by the L 2 norm, over all grid points, of the difference between the computed and 
“exact” solution at some final time. The final time corresponds to the movement of the shockwave 
approximately 100 shock thicknesses. A quasi exact solution is found bv running ARK4(3)6L[2]SA 
at an order of magnitude finer time step than any used in the grid refinement study. At e = 10 -6 , 
error increases but the order of accuracy remains nearly three. Two seperate measures of efficiency 
may be considered: accuracy and stability efficiency. Accuracy efficiency determines the work required 
to obtain some chosen error tolerance. We define work as the number of implicit solves required for 
the integration without regard to Newton iteration count. O 11 error versus work plots, the five stage 
methods of Ascher (4, 4, 3) and Fritzen FW53 are least efficient. LIRK3 is the least accurate of the four 
remaining methods on this particular problem, followed by Ascher (3, 4, 3). The most accuracy efficient 
methods are ARK3(2)4L[2]SA and Ascher (2,3,3), as shown in figure 2. There was no evidence of the 
coupling stability term causing a problem with Ascher (2,3,3). When accuracy is sufficient, one simply 



Figure 2: Error versus work for AR.K3(2)4L[2]SA and Ascher (2,3,3) in the presence and absence of 
stiffness. 

seeks the largest stable time step. The limiting time step might be due to ERK linear stability boundary 
or iterative convergence problems with the Newton’s method. In the absence of any stiffness leakage 
or convergence difficulties, one may compute a theoretical stability based efficiency of the methods by 
considering the inviscid and viscous CFL numbers normalized by the number of implicit stages. This 
minimum work point would correspond to the upper right limit of the lines on an error versus work 
plot. Table 14 compares the maximum time step per unit work for the first- and second-order schemes 
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of Ascher et al.,° six third-order schemes, two fourth-order schemes, and the one fifth-order method. 
ARK5(4)8L[2]SA has marginal stability efficiency for both real and imaginary lionstiff eigenvalues. Both 
LIRK3 and LIRK4 may not be efficient for relatively stiffer convective eigenvalues. In compressible 


Method 

CFLi nv iscid/ Sj 

CFL viscous / Sj 

Ascher (1, 2, 1) 

0.435 

0.315 

Ascher (2, 3, 2) 

0.500 

0.250 

Ascher (2, 3, 3) 

0.435 

0.315 

Ascher (3,4,3) 

0.473 

0.233 

Ascher (4,4, 3) 

0.195 

0.135 

Calvo LIRK3 

0.023 

0.183 

Fritzen FW53 

0.218 

0.158 

ARK3(2)4L[2]SA 

0.413 

0.307 

Calvo LIRK4 

0.032 

0.176 

ARK4(3)6L[2]SA 

0.402 

0.212 

ARK5(4)8L[2]SA 

0.061 

0.096 


Table 14: Idealized stability based efficiencies of first- through fifth-order methods in the absence of 
stiffness leakage. 


flows, it is likely that the inviscid limit is more relevant to observed stability efficiency. Which of these 
two limits is more important during stiffness leakage is not clear. Comparing tables 11 and 14 at 
e = 10° suggests that the when stiff real eigenvalues leak, the real axis (viscous) stability efficiency may 
be more important. Using equilibrium initial conditions and exothermic chemistry, Ascher (2,3,3) was 
found to be more stability efficient than ARK3(2)4L[2]SA at e = 10 -6 ; however, ARK3(2)4L[2]SA may 
possess some intrinsic efficiency advantage over Ascher (2,3,3) in having a smaller value of 7: 0.4359 
versus 0.7887. Modified Newton iteration would presumably converge faster. 

At fourth-order and above, for DIRK-based IMEX methods, we are only aware of the LIRK4 method 
attributed to Calvo et al. 10 and the methods that we have generated. Calvo LIRK4 was constructed 
by adding an ERK to the zero-padded SDIR.K4. 2 ' Both the implicit and additive methods fully damp 
stiff scaled eigenvalues. Large leading-order error of the ER.K dominates the leading-order error of the 
IMEX method, as shown in appendix B. The method may conserve quadratic first integrals poorly and 
has a relatively small convective stability limit for the ER.K. ARK4(3)6L[2]SA results from an extensive 
examination of possible approaches to fourth-order methods and is described in section 5. 

Two particular methods used (7(2,?)^, (7(2. *)W, and 7 = 1/4. One had A^ 5 ) = 0.01224 and the 
other Al ' = 0.01284, but used simplifying assumption (7(3, iy- 1 * also. With nearly identical overall 
error and slightly worse internal stability values, both methods behaved nearly identically. The next 
test was to compare use of simplifying assumption <77(2, i)^ with no assumption for the ERK. Leading 
order error for the methods were Al 5 ^ = 0.01224 and A = 0.01461, respectively. Without (7(2,*)^, 
the method becomes progressively less accurate relative to the other method as stiffness increased. To 
test the effect of coupling stabilty, two methods were constructed employing (7(2,*)^, (7(2, i)^, and 
7 = 1/4. Each satisfied £>60 = Pso = p . 51 = 0, but the second additionally satisfied £>41 = £>42 = 0. The 
second method has a much smaller linear stability domain for its ERK and a principal error norm of 
Al 5 ) = 0.03542, three times that of the first. It appeared to be slightly less prone to stiffness leakage 
than the first scheme. Possibly the more highly optimized linear stability domain of the first method 
was more susceptable to degradation. Interestingly, in the presence of stiffness, several fourth-order 
methods retain a smooth reduction in error as work is increased while others do not. Methods using 
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C( 2, 2')^ appear smooth while those using no simplifying assumption on their explicit method appear 
rather jagged, e.g., LIRK4 and the five-stage, 4(2) pair that was constructed. One may also compare the 
degree of error increase with the addition of stiffness. LIRK4 has A( s ) = 0.03919 and uses no simplifying 
assumptions whereas one of the test methods has A* 5 ) = 0.03542 but makes use of (7(2, *) [B1 , C(2,*) [i] . 
Their nonstiff accuracy efficiencies are quite similar but in the presence of stiffness, LIR.K4 shows not 
only a more dramatic increase in error but more order reduction. Since LIRK4 does not use C( 2, 
and C(2,i)W, it is not unreasonable to attribute this difference to the stage-order of the implicit method. 
Finally, results of various methods on this test problem often correlated with the leading-order error 
term of the entire method, Figure 3 compares Calvo LIR.K4 and AR.K4(3)6L[2]SA at stiffness 



Figure 3: Error versus work for ARK4(3)6L[2]SA and Calvo LIRK4 in the presence and absence of 
stiffness. 

extremes showing that ARK4(3)6L[2]SA is not only more accurate but increasingly so as the stiffness 
is increased. 

The only fifth-order method that we are aware of, AR.K5(4)8L[2]SA, offered few design options. 
Assumptions C( 2, i)^ and (7(2, «)^ were essential. To reduce embedded order conditions, one column 
each of both Butcher arrays was replaced with zeros. Doing this made it only reasonable to include 
one (7(3, i) assumption; hence, (7(3, i)^ was used along with D(l,j)^ to reduce the number and 
complexity of the remaining order conditions. Selection of 7 was problematic. Schemes were designed 
using 7 = 0.145, 0.1725, 0.185, and 0.205. Larger values appeared to work better. The resulting 
method, ARK5(4)8L[2]SA, has a relatively small linear stability region for its ERIv. In the absence of 
stiffness it is found to be fifth-order, but in the presence of strong stiffness it order-reduces in the same 
manner as ARK4(3)6L[2]SA. Given the behavior of ARK5(4)8L[2]SA on these tests, it is probably the 
best choice for situations where both mild stiffness and tight error tolerances are present. Otherwise, 
ARK4(3)6L[2]SA might best be employed. Finally, comparing efficiency of methods of all orders at 
£ = 10° in figure 4 and at £ = 10 -6 in figure 5, one may conclude that first-order methods are ill- 
advised. Ascher (2,3,2) can provide a much more accurate solution for an identical cost. At any level of 
stiffness, Ascher (2,3,2) is very stability efficient. At tighter error tolerances, ARK4(3)6L[2]SA would 
appear to be the most efficient high-order IMEX additive Runge-Kutta method of w r hich we are aware 
for this problem. 
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Figure 4: Error versus work for first- through 
fifth-order methods in the nonstiff case. 



Figure 5: Error versus work for first- through 
fifth-order methods in the stiff case. 


A limited attempt is made here to quantify the behavior and efficiency of the multistep SBDF 
schemes as they are the chief alternative to the IMEX Runge-Kutta schemes. Multistep schemes 
solve only one nonlinear set of equations per time-step. As a result, they can potentially achieve 
great efficiency compared with multistep Runge-Kutta formulations. Like most ARK 2 methods in 
this study, the SBDF methods lack error control. Additionally, SBDF methods are not self starting, 
require fixed time steps, and the implicit formulations do not possess A-stability beyond second order 
accuracy. Second-, third- and fourth-order SBDF schemes of Ascher et al. 5 achieve design accuracy on 
the reacting shock wave problem independent of the stiffness level, although the leading order constant 
appears to depend on stiffness. One may compare the various Runge-Kutta and SBDF methods in 
terms of either accuracy or stability efficiency. In terms of accuracy, Ascher (2,3,2) is substantially 
more efficient than SBDF2; however, at higher-order, the SBDF methods are more efficient on this 
problem. Figures 6 and 7 compare third- and fourth-order multistep and Runge-Kutta methods with 
and without stiffness. It should be remembered that chemical reaction rate terms generally give rise 
to real eigenvalues. The maximum time step for the SBDF formulations is strongly dependent on the 
location of the scaled stiff eigenvalues, and in particular, whether they fall in the unstable lobes of 
the implicit BDF.3 and BDF4 operators. Scaled eigenvalues on or near the negative real axis are well 
suited for implicit BDF operators, while eigenvalues near the imaginary axis are not. Conversely, RK 
schemes degrade in accuracy with increased stiffness due to their lower stage-order. 

8.3 Order- Reduction 

Each of the four test problems in this paper are examples of singular perturbation problems. 25, 2 ' 
They are ODEs characterized by a stiffness parameter, e. As e decreases, the ODE problems gradually 
transition in behavior toward index- 1 DAEs. For Runge-Kutta methods, accompanying this transition 
is an order-reduction phenomena where the observed convergence rates of the methods fall below the 
classical order of accuracy. Some differential variables transition to algebraic variables, displaying 
different convergence rates. Hairer et al. 25 ’ 2 ' determine the convergence rates of SDIRK methods with 
and without the stiffly accurate assumption. Global error for both differential and algebraic variables 
are of the form e g i 0 b a i = ci(A t) a + C 2 e(At) 13 for e < Const. (At). Independent of stiff accuracy, SDIRK 
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Figure 6: Efficiency comparison of third- and Figure 7: Efficiency comparison of third- and 
fourth-order methods with e = f0°. fourth-order methods with e = 10~ 6 . 

methods have a = q and /3 = q so + 1 for the differential variable, where q and q so are the classical and 
stage-orders of the method. In the case of algebraic variables, a = q so + 1 and C 2 = 0 for non-stiffly 
accurate methods but a = q and /3 = q so in the stiffly accurate case. Practical experience shows that 
the order reduction is problem dependent and may not be as severe as the theoretical estimate. 

Although theoretical bounds exist for many implicit Range- Kutta methods applied to singular 
perturbation problems, little exists for IMEX methods.' 3 No attempt is made to theoretically predict 
the form of the global error of IMEX AR.K 2 methods, rather we shall use the values and form articulated 
above for SDIRK methods as guidance in empirically estimating the values of the respective exponents. 
The ultimate goal of this order-reduction study is to establish its severity for CDR problems using stiff 
chemical kinetic mechanisms. Both the vdP and CDR problems will be used to this end. A cursory 
examination of the findings from the previous section might lead one to the incorrect conclusion that 
no order-reduction exists in the CDR problem. Order reduction is easily identified in simple model 
problems. Thus, we begin our study with three stiff singular-perturbation model problems. We establish 
the accuracy of the new methods on these problems and compare them to existing SDIRK schemes in 
the literature. Attention is then focused on the reacting shock wave problem where its order-reduction 
characteristics are demonstrated. 

8.3a Order- Reduction on Model Problems 

All previously mentioned numerical schemes were run on Kap’s problem, van der Pol’s equation, 
and Pareschi and Russo’s problem. In each case, fully implicit and IMEX formulations were compared 
to assess the effects of partitioning. Order reduction is observed for all ARK 2 schemes whose classical 
order is greater than two, but is not observed for the SBDF formulations. The general nature of 
the order reduction is similar for all three problems although the degree of reduction varies between 
problems. Since van der Pol’s equation exhibits the greatest order-reduction, it is chosen as the testbed 
to compare the accuracy of all schemes. The time interval chosen for these studies is 0 < t < 0.5. 
Figures 8 and 9 show representative results of a temporal refinement study at various levels of the 
stiffness parameter e for van der Pol’s equation obtained using ARK4(3)6L[2]SA. From equation (61), 
y 1 is the differential variable, while y? transitions from a differential variable to an algebraic variable 
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Figure 8: Nonstiff (differential) error in the Figure 9: Stiff (algebraic) error in the van 

van der Pol equation as calculated with the der Pol equation as calculated with the 

ARK4(3)6L[2]SA scheme. ARK4(3)6L[2]SA scheme. 


as the stiffness is increased. The convergence rate of the differential variable is nearly fourth-order, its 
classical order, and is reasonably smooth. Convergence behavior of the algebraic variable, however, is 
jagged and departs significantly from the expected design accuracy. This departure from the design 
accuracy occurs most dramatically for intermediate values of the stiffness parameter e. 

In figures 10, 11, and 12, the convergence rates of the ARK3(2)4L[2]SA, ARK4(3)6L[2]SA, and 
ARK5(4)8L[2]SA methods on van der Pol's equation are plotted versus the stiffness parameter e. 
Convergence rates are calculated by a least-squares fit of the data, at each e presented in figure 8. This 
procedure is at some variance with the fact that error versus stepsize lines are generally composed of 
two lines of differing slopes. It is adopted, nonetheless, because the dual slope lines are not discernable 
from the jagged data. ESDIRK method results, for which theoretical estimates of order reduction exist, 
are included for comparison purposes along with the IMEX values. Certain general trends appeared 
across each test problem. At values of £ « 10°, the observed convergence rate was equal to the classical 
order of each method. 

All curves show order reduction to some degree for intermediate values of the parameter s. The 
£ — * 0 limit is again characterized by a uniform convergence rate, although it is not generally the design 
rate of the method. Both formulations, the ESDIRK and the IMEX, order-reduce for the algebraic 
variables considerably more than for the differential variables. The general trends presented in figures 
10, 11, and 12 are summarized in table 15 along with the Ascher (2,3,3) scheme, Fritzen FW53, and 
Calvo LIRK3 and LIRK4. Based on results from the three model test problems, we determine the 
leading order truncation terms for each method for the cases where £ < Const. (At). In practice, this 
is an awkward and inexact procedure due to the extreme jaggedness of the convergence plots. With 
this in mind, some methods showed problem dependent convergence behavior. The exact nature of a 
method’s order reduction depends on the relative size of the parameters £ and At. The £ <C At limit 
yields the £(At)^ contribution to the convergence rate for the method (see figure 10 where £ = 10~ 8 ). 
The £ At limit gives the classical order for all methods as £ is not a small parameter in this case. 

Onset of order reduction is observed for the cases where £ ps At. By fixing £ and varying At in the 


32 






Figure 10: Convergence rates 
of ARK3(2)4L[2]SA on the 
vdP equation as a function of 



Figure 11: Convergence rates 
of ARK4(3)6L[2]SA on the 
vdP equation as a function of 



Figure 12: Convergence rates 
of ARK5(4)8L[2]SA on the 
vdP equation as a function of 
£. 


Method 

ESDIRK 

Differential 

ESDIRK 

Algebraic 

IMEX 

Differential 

IMEX 

Algebraic 

Ascher (2,3,3) 
Cal vo LIRK3 
Fritzen 

ARK3(2)4L[2]SA 

(At) 3 + £(At) 2 
(At) 3 + £(At) 2 
(At) 3 + f( At ) 2 
(At) 3 + f(At) 3 

(At) 2 + s(At) 1 
(At) 3 + e( At) 1 
(At) 3 + e(At) 1 
(At) 3 + e(At) 2 

(At) 3 + s( At) 2 
(At) 3 + e(At) 2 
(At) 3 + £(A t) 2 
(At) 3 + £(At) 2 

(At) 2 + £( At) 1 
(At) 2 + £( At) 1 
(At) 3 + £( At) 1 
(At) 2 + £( At) 1 

Cal vo LIRK4 
ARK4(3)6L[2]SA 

(At) 4 + f(At) 2 
(At) 4 + £( At) 3 

(At) 4 + s(At) 1 
(At) 4 + f( At) 2 

(At) 4 + £(At) 2 
(At) 4 + s(At) 2 

(At) 2 + £( At) 1 
(At) 3 + £( At) 1 

ARK5(4)8L[2]SA 

(At) 4 + e(At) 3 

(At) 4 + e( At) 2 

(At) 4 + s(At) 2 

(At) 3 + £( At j 1 


Table 15: Estimated convergence rates of differential and algebraic variables in e g i 0 b a i = Ci(At)“ + 
C 2 £{&t)P form for several ESDIRK and IMEX ARK 2 methods for e < Const. (At). 
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£ At limit, a careful study of order- reduction may be conducted. Typically the convergence rate for 
these studies has a slope discontinuity, that can be used to identify the leading order error terms. For 
the IMEX operators, convergence behavior is erratic in this limit, often not monitonically decreasing 
with time step. Convergence behavior for the third and fourth-order implicit schemes, presented in 
table 15, agrees with the theoretical estimates for SDIRK schemes given above. Incorporating a stage- 
order of two in the implicit scheme does not translate into higher order-of-accuracy for the differential 
variables for any of the IMEX schemes, but it increases the convergence rate for the ESDIRK alone. 
Also, the increased stage order increases the accuracy of the algebraic variable in the fourth- and fifth- 
order cases. The 4tli-order asymptotic convergence of scheme ARK5(4)8L[2]SA and its ESDIRK makes 
the overall behavior of the fifth-order scheme very similar to that of the 4th-order scheme. We offer no 
explanation for this behavior. 

8.3b Order- Reduction on CDR problems 

We now extend the previous study on order reduction to CDR problems using the reaction inducing, 
propagating shock wave problem. Computations use exothermic chemistry and equilibrium initial 
conditions. As the reaction rate is introducing the stiffness and only the species continuity equations 
explicitly contain the reaction rate, we focus on these equations. From Eq. (63), the species equations 
are given by 


d(pYj) 

dt 


d(puYj) 

dx 


„ d 2 Y t . 
+ PDi ~d^ + UJi 


Ans.e' T dy , 


(70) 


where Y) = \Y}{ 2 ,Yo 2 ,Yy{ 2 OiYoy{,Yn 2 } and F nbt denotes the sum of convective and diffusive terms for 
species i. Because of overall continuity, it is not necessary to solve a differential equation for Yn 2 • 
Nonstiff convective and diffusive terms are integrated with the ERK while all reaction rate terms are 
integrated by the ESDIRK. Using concentrations, C l = pY t /W t , and our simplified reaction mechanism, 
the full species equations appear as 


pYn 2 


-^ns,H 2 


' IEH 2 aexp(-^/^ [k(-Cu 2 Co 2 + Cg H ) - Ch 2 C§ h ] 
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where e = k^/k^ and a is some real constant. Fast reactions involve the term 7 and are present in 
the reaction terms for the species H 2 , O 2 and OH. Variable O 2 is purely a fast reaction. The other 
two variables, H 2 and OH, involve both fast and slow reactions, while the variable H 2 0 is entirely 
a slow variable. In general, it is difficult to identify slow (differential) and fast (algebraic) variables 
and hence treat each appropriately. For this simplified reaction system, however, the differential and 
algebraic variables can readily be identified. Defining the new variables £1 = Wo 2 Yh 2 - Wh 2 Vq 2 and 
£ 2 = ll o 2 Voh + 2 !I ohVo 2 , the system may now be written as 



1 
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Wo 2 F ns ^2 ~ WH 2 ^ns,0 2 

d 

pYo 2 


-^ns,0 2 

dt 

pYn 2 o 


Fjxs,H2 


Pi 2 


_ VFq 2 -^ns,OH - 2l4 7 OH-f 1 ns,0 2 


Wu 2 Wo 2 ol exp (~ T o/ T f [ — C'h 2 Cq H ] 

Wo 2 ae^~To/Tf [1(-C H2 Cc 2 + C 2 oh )} 

2W H2 oo:ex p YTo/Tf [6'h 2 6'o II ] 

. 2WouWo 2 aex V (~ T ^ T f [-Cn 2 Cl n ] . 


With this new system, it is more clearly seen that at high levels of stiffness, Yo 2 will be an algebraic 
variable while £| , Vh 2 o, £25 and V n 2 } are the differential variables. 
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Figure 13: Convergence rates of differential and algebraic variables on CDR problem using the new 
ARK 2 methods. 

Figure 13 compares the convergence behavior of the three new IMEX ARK 2 schemes on the reacting 
shock wave problem. A representative differential variable, H 2 O, and the algebraic variable, O 2 , are 
presented for ARK3(2)4L[2]SA, ARK4(3)6L[2]SA, and ARK5(4)8L[2]SA. Convergence rates are again 
determined from a least-squares fit of a convergence study at each value of the parameter s. Interme- 
diate values of the parameter s, again, produce the most order-reduction in both the differential and 
algebraic variables. 

The algebraic variable order-reduction in this CDR. problem is remarkably similar to that observed 
in all three singular perturbation model problems. This degradation in accuracy for the algebraic 
variable, however, does not dramatically degrade the overall accuracy of the C'DR problem when error 
is considered as the L 2 norm of the difference between the computed and exact solutions over all grid 
points. Temperature, which is a combination of differential and algebraic variables, converged at a rate 
slightly lower than the differential variables. As stiff modes are unnecessary to resolve for accuracy 
purposes and algebraic variables arise from high stiffness, it may not be surprising that the lower 
convergence rates of algebraic variables appear to weakly affect temporal error. Another explanation 
for the benign role of order-reduction in the present CDR problem is that the one algebraic variable 
is only weakly coupled to the rest of the system and does not greatly influence the solution accuracy 
of the other six variables. It is not clear if this may be generalized to all or most reacting flows, 
however, AR.K 2 schemes are likely to experience significant order-reduction on problems where the 
algebraic component of the error plays a dominant role. In this scenerio, the SBDF schemes, which do 
not experience order-reduction, are likely to have a clear efficiency advantage over the IMEX Runge- 
Ivutta schemes provided the stiff eigenvalues are predominately real. A case by case study is probably 
necessary to definitively answer whether order- reduction is an important issue. 

8.4 Error Control 

Choosing a practical error controller for the current IMEX methods is problematic. Advanced 
controllers designed for explicit and implicit methods are constructed based on different criteria. IMEX 
schemes, being combinations of each, represent a new challenge for error controllers. Beyond this, 
with increasing stiffness, controllers additionally confront order- reduction as well as emerging algebraic 
variables. With this in mind, we test four general appoaches: the I , PI-, PC, and P ID -controllers. 63 
The 1-controller is appropriate for either implicit or explicit methods. PI- and PID-controllers are 
advances over I-controllers for explicit methods. PC-cont rollers have been designed for the unique 
dynamics of an implicit method. 
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Van der Pol’s equation provides a challenging test for the error control capabilities of the new 
additive schemes. Over one period, the vdP solution has two temporal boundary layers having thickness 
related to e -1 . To maintain accuracy, the embedded method and controller must sense the layers and 
adjust the time step accordingly. Significant variations were observed between different controllers. For 
example, the PC-controller is ineffective. It produced a strong step-change instability characterized 
by large time-step changes in portions of the temporal cycle where no adjustments were necessary. 
Controllers designed for implicit methods appear inappropriate. Performance of the simple I-controller 
is better, yet marginal. Both the PI- and PID-controllers are able to guide the integration through the 
temporal boundary layer with reasonable efficiency 7 . With better SC-stability properties and similar 
characteristic roots, the PID-controller behaved best. As the PI- and PID-controllers are designed for 
the dynamics of explicit methods, it appears that controlling the IMEX method on these problems 
is largely a function of controlling the explicit method. Further, it appears that the behavior of the 
controller at the stability boundary is most important. That these controllers worked well is surprising 
considering that in the presence of large stiffness, the algebraic variables dominate solution accuracy in 
the three model problems and order reduction is present. The generality of these findings is not clear. 

In a second test of the error controllers, the propagating reacting shock wave problem is computed. 
Testing methods on this problem is rather difficult because the flow contains no transients, yet if one 
specifies an exothermic reaction system, a nonequilibrium initial condition, and e = 10 -6 , a highly 
transient problem ensues. This same problem is severe enough to break all ARK 2 schemes when used 
in fixed stepsize mode and large stiffness. Conclusions are similar to those drawn from van der Pol’s 
equation. The only controllers capable of guiding the integration out of the flow equilibration phase 
at all stiffnesses are the PI- and PID-controllers. At low and high stiffness, requested and resultant 
error are well correlated. Stiffness affects the relation of predicted and actual error but the controller 
remains useful and is remarkably insensitive to the stiffness even in the case of order reduction. It 
is less surprising that the explicit-based controllers perform adequately on the CDR problem as the 
algebraic variables are of secondary importance to solution accuracy. 

We do not offer any theoretical explanation why the PI- and PID-controllers work fairly well on 
these problems. To maintain constant controller gain during order-reduction, p should presumably be 
reduced. It is not reduced in these tests. Perhaps the essential feature of controlling IMEX methods is 
coping with scaled eigenvalues at the stability boundary of the explicit method, a task best suited to 
the PID-controller. 

8.5 Dense Output 

The dense output for the three new schemes is tested on the reacting shockwave problem. An 
equilibrated solution is established at a time, t = at £ = 10°, and is used as the initial condition for 
the study. The initial condition is then advanced one time step to fill all function registers. Interpolation 
and extrapolation are done at points preceding and following f re f + At. The dense output is then 
compared to an “exact” solution obtained with a separate run beginning with the initial condition 
using At/10 stepsizes and run to the dense output times. A refinement study is performed using one 
timestep in the variable At to determine the local order of accuracy of the dense output. Note that 
the nature of the refinement study in the variable At returns the local error of the dense output or 
the global error plus one. Table 16 summarizes the observed local errors, S ^ t , and convergence rates 
from a study using the third-order (p* = 3) formula associated with the ARK4(3)6L[2]SA scheme. The 
interpolated and extrapolation values are at |(At) and |(A t), respectively. 

Design order is asymptotically achieved in both modes. Note that the extrapolated data are one and 
one half orders less accurate than the interpolated data although their respective orders-of accuracy 
are similar. The efficacy of extrapolation decays rapidly with distance. Similar results showing design 
order dense output were obtained for the ARK-3(2)4L[2]SA and ARK5(4)8L[2]SA schemes. A final test 
of the dense output was performed on both van der Pol’s and the CDR equations. Extrapolation mode 
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was used to predict the starting values of the newton iteration from data at the previous timestep. 
A uniform speedup of the iteration was observed at all levels of stiffness, e = 10° through s = 10 -6 , 
indicating the efficacy of the extrapolation. 


Table 16: Convergence rate and local error of the p* = 3 dense output for interpolation and extrapo- 
lation as calculated with the ARK4(3)6L[2]SA scheme. 


At 

^(Ai),Int 

Order 

At), Ext 

Order 

0.9 

-5.218 


-3.599 


0.6 

-5.790 

3.241 

-4.229 

3.579 

0.5 

-6.070 

3.541 

-4.525 

3.745 

0.4 

-6.426 

3.679 

-4.896 

3.823 

0.3 

-6.900 

3.791 

-5.382 

3.891 

0.2 

-7.583 

3.879 

-6.077 

3.945 

0.1 

-8.770 

3.944 

-7.275 

3.981 


9. Conclusions 

Additive Runge-Kutta (ARK) methods are investigated for application to the spatially discretized 
one- dimensional convection-diffusion-reaction (C'DR) equations. First, accuracy, stability, conserva- 
tion, and dense-output are considered for the general case when N different Runge-Kutta methods 
are grouped into a single composite method. Comparing the N = 3 and N = 2 cases for CDR ap- 
plications, N = 2 methods are chosen. Then, implicit-explicit, N = 2, additive Runge-Kutta (ARK 2 ) 
methods from third- to fifth-order are presented. Each allows for integration of stiff reactive terms 
by an L-stable, stiffly-accurate ESDIRK method while the nonstiff convection and diffusion terms are 
integrated with a traditional ERK method. Coupling error terms are minimized by selecting identical 
abscissae and scheme weights for each method and are of equal order to those of the elemental methods. 
Both ARK 2 and ESDIRK methods have vanishing stability functions for very large values of the stiff 
scaled eigenvalue, z^ — »• — 00 , and retain high stability efficiency in the absence of stiffness, z^ — > 0. 
Extrapolation-type stage- value predictors are provided based on dense-output formulae. Dense output 
stability functions have minimized values for 0 > 1 and — >■ — 00 . Optimized methods minimize 

both leading order ARK 2 error terms and Butcher coefficient magnitudes as well as maximize con- 
servation properties. Numerical tests of the new schemes on a CDR problem show negligible stiffness 
leakage and near classical order convergence rates. Third- and fourth-order SBDF methods are slightly 
more efficient than the IMEX ARK 2 schemes but do not include error estimation and stepsize control. 
Tests on three simple singular perturbation problems reveal similar and predictable order reduction 
for the Runge-Kutta. methods but no order reduction for the SBDF methods. Order reduction of 
ARK 2 schemes is worst at intermediate stiffness levels. Estimated convergence rates for differential 
and algebraic variables generally coincide with that predicted by theory. A reinspection of differential 
and algebraic variables on the CDR problem shows similar behavior. Error control is best managed 
with a PID-controller, indicating that ERK stability is the overriding issue in controlling error. Dense 
output is useful both in interpolation and extrapolation. While results for the fifth-order method are 
disappointing, both the new third- and fourth-order methods are at least as efficient as existing ARK 2 
methods while offering error control and stage- value predictors. 
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Appendix A - Runge-Kutta Order Conditions 

Equations of conditions for 1-trees up to sixth-order accuracy are given by 
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General equations of condition for 2-trees up to fifth-order for one root node type are provided below. 
Coefficients of the two methods are distinguished by case; (a 8 -j, c;) and ( Aij,B{ , C,-). Only half of the 

actual order conditions are given. The other half may be obtained by taking each condition below and 
replacing each b{ with i?;. 
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Appendix B - ERK/IRK Additive Runge-Kutta Methods 
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Appendix C - ERK/IRK Additive Runge-Kutta Embedded 
Methods 
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Appendix D - Additive Runge-Kutta Scheme Coefficients 
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ARK3(2)4L[2]SA - Second-Order Dense Output 
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